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Abstract 



The problem of modeling forest tree growth curves with an artificial neural network (NN) 
is examined. The NN parametric form is shown to be a suitable model if each forest tree 
plot is assumed to consist of several differently growing sub-plots. The predictive Bayesian 
approach is used in estimating the NN output. 

Data from the correlated curve trend (CCT) experiments are used. The NN predictions 
are compared with those of one of the best parametric solutions, the Schnute model. Analysis 
of variance (ANOVA) methods are used to evaluate whether any observed differences are 
statistically significant. From a Frequentist perspective the differences between the Schnute 
and NN approach are found not to be significant. However, a Bayesian ANOVA indicates 
that there is a 93% probability of the NN approach producing better predictions on average. 
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Chapter 1 



Introduction 



Growth curves are found in many different areas of study, e.g. in animals, plants, bacteria, 
fishes and populations. Their study is important in understanding what they are affected 
by and in predicting future values. 

In the field of forestry much effort has been expended in finding mathematical models 
to describe the growth of trees (Bredenkamp and Gregoire 1988, Beber and Wild 198S|, 



Falkenhagen 1997). Models such as the logistic function have been proposed for predicting 



the average tree diameter at breast height (DBH) for a plot of forest trees: 



y(i) = 



1 — exp(— ut + a] 
logistic(i; u, w, a) 



where y{t) is the average DBH at age t and u, v and a are the model parameters. For the 
model to be physically realistic: 



(1.1) 



t > 0. 



An advantage of such parametric models is that their parameters are easy to interpret, e.g. 
V will be the maximum average DBH attainable. The disadvantage of such models is that 
they are only appropriate for modeling a very limited family of input/output mappings. 



Multilayer perceptron artificial neural networks (NNs) ( Haykin 1994 ) provide a more 
flexible method of nonlinear regression. A general functional form of a NN for a one dimen- 
sional input to output mapping is given by: 



(1.2) 



Nh 



y{t) = ^ logistic(t; Ui, Vi,a.i) + b. 



The larger Nh is chosen to be, the more flexible the model. Hornik et al. (199C ) have 



shown that equation (1.2) has fairly general function approximation qualities. In the case 
of modeling average forest tree growth, the NN model has a particularly suitable form. If 
the forest tree plot can be assumed to consist of Nh groups of differently growing trees and 
each group's average is modeled using the logistic function, then the NN functional form 



1 



follows. However, the NN model does not usually contain the physical constraints mentioned 



in equation (1.1) 



Thus, the NN model provides for heterogeneity in the growth of a forest tree plot. 
Unlike most NN applications, in mean forest tree growth modeling the NN functional form 
has some justification. 

1.1 Artificial Neural Network Parameter Estimation 

In order for the NN to have enough flexibility to fit a wide range of growth curves, Nh has 
to be made fairly large. But a large Nh implies many model parameters. The more model 
parameters, the more sensitive the solution is to any statistical variability in the data. This 
is known as the bias/variance dilemma. 

Bayesian estimation provides a way of achieving a low bias without paying the price of a 
high variance. Predictions of unmeasured growth values are made by taking a weighted 'sum' 
of the predictions provided by all possible values of the parameters. Given N input/output 
pairs, V = {ti,yi; . . . ;tjv, J/at}, a new measurement, ijn+i, is estimated by: 



yN+i= / y{tN+i;e)p{6\V)de 

where 9 is the set of NN parameters and Vg is their domain. The weight of each prediction 
is given by the probability density function (pdf ) of the network parameters given the data, 
p(9\'D). This pdf can factorized as follows: 

p{e\v) (X p{v\e)p{e) 

where p{9\'D) is known as the likelihood and expresses how the pdf is affected by the available 
data. The p{9) component is known as the prior and expresses data (I?) independent 
knowledge about the model parameters. The prior can be used to include appropriate 
smoothness constraints which reduce the variance of the estimate. 



1.2 Model Evaluation 

When deciding which model is better at describing the process that generated a particular 
set of data, it is preferable to test the model on different data than it was trained on. 
Otherwise, the performance of each model is likely to be optimistically biased. 

Using only one training set and test set to compare regression methods can be mis- 
leading. How well each method does will depend on the particular training and test cases 
used. The methods should be compared using many different training sets and test cases. 
Statistical hypothesis testing can then be used to determine which method is on average the 
best. 

1.2.1 Correlated Curve Trend Experiments 

A set of forest tree growth data suitable for the comparison of regression methods can 
be obtained from the results of what are known as the "correlated curve trend" (CCT) 



experiments ( O'Connor 1935| ) . The growth of several different plots of trees with different 



initial and growing conditions was monitored. The growth measurements for the different 
plots can be used to provide an estimate of how well the NN approach performs on average 
in comparison with a parametric regression approach. 

1.3 Objectives 

In this report a survey of the available methods of forest tree growth curve modeling and 
Bayesian artificial neural networks (BNN) will be given. Statistical hypothesis testing will 
be used to compare the BNN approach with standard parametric models on the forest tree 
growth modeling problem. 

Another of the aims of this research report is to evaluate the BNN approach on a 
practical, real world problem where a relatively small amount of data is available. 

1.4 Outline 

In Chapter the previous literature on forest tree growth curve modeling is reviewed. The 
CCT experimental data are examined and the objectives of the curve fitting procedure are 
defined. The Schnute solution proposed by Falkenhagcn (1997| ) is also discussed. 



In Chapter y the Bayesian methodology used in the rest of the report is reviewed. 
Many of the relevant results are derived from first principles. The problems of defining 
prior distributions and the derivation of predictive distributions are discussed. 

Chapter surveys the relevant NN literature. The bias/variance dilemma is explained. 
The Bayesian treatment of Neal (1996|) is summarized. A new justification for the prior 



distribution assigned to the network weights is given. 

Chapter || reviews an analysis of variance (ANOVA) scheme, introduced by Rasmussen 



(1996 ), for comparing regression methods. A Bayesian hierarchical solution is also discussed. 
In Chapter || the Bayesian neural network (BNN) scheme is applied to extrapolating 



the CCT data. The Frequentist ANOVA approach of Rasmussen (1996 ) and a hierarchical 
Bayesian ANOVA are used to compare the statistical significance of the Schnute and BNN 
results. 

An overview of the results obtained in this research report is presented in Chapter m. 
The scope and limitations of the results are discussed. 



Chapter 2 

Forest Tree Growth Modeling 



If physics has its laws, biology has its variety. - G. A. Dover. 



2.1 Background 

There is a long history of forest tree growth modeUng, from the first yield tables published 



over 200 years ago, to the recent Bayesian treatments of growth and yield models (Vanclay 



1995, Green and Strawderman 1996). Models help in forecasting timber yields, identifying 



appropriate treatments, planning how densely to plant trees together, deciding when to 
harvest and in monitoring the current state of a forest. They also help in determining the 
sustainability of various silviculture practices. 

Vanclay (1995) has given a synthesis of the models and methods for tropical forests. 
An important aspect of tropical forest modeling is whether the timber harvesting is sustain- 



able (Vanclay 1994). Oshu (1991) uses a matrix model to predict long term tropical rain 
forest growth, in which matrix eigenvalues are used to estimate the intrinsic rate of natural 
increase. Methods of assessing the usefulness of permanent sample site databases are given 



Vanclay et al. (1995]) 



Bayesian techniques have been used to estimate the parameters of a growth and yield 
model for slash pine plantations (Green and Strawderman 1996). Posterior probability 
distributions were found for parameters such as number, volume and diameter of plantation 
trees. Zellner's (199q) Bayesian method of moments was used to avoid having to make 
any assumptions about the form of the likelihood function. Another Bayesian paper is 



Green et al. (1994 ) where Bayesian estimation is used to fit the three parameter WeibuU 
distribution to some tree diameter data. It is shown that the Bayesian solution avoids the 
negative location parameter estimates which plague the maximum likelihood solutions. 

In a series of articles, Guan and Gertner (1991a q;n; 1995) used an artificial neural 
network (see Chapter ^) to model forest tree mortality in terms of diameter at breast height 
(DBH) and the annual increment in DBH. 



2.2 Correlated Curve Trend Experiments 

The growth of a tree can be afFected by competition from neighbouring trees for the available 
resources of sunlight, moisture, root space and soil nutrients ( Vanclay 1995| ). The degree of 
crowding has a considerable effect on the mean tree size. 

O'Connor (1935) has examined the question of how the crowding of trees effects their 
growth. There are two components to this problem: 

1. How closely the trees are planted together, known as the espacement. 

2. What thinning^ strategies are employed. 

There are a number of different qualities that a forester might consider when determining 
the optimum strategy: 

1. The total volume of production, e.g. for pulp production. 

2. How quickly the trees will grow. 

3. The distribution of tree sizes. 

The effects of different thinning regimes can be determined by keeping all other relevant 
factors constant and just varying the thinning regime. Generally the main factors in de- 
termining tree growth are the species of tree and the location or site where the trees are 
growing. Thus to compare different thinning regimes, the same forest tree species is planted 
on a site which is as uniform as possible. 



The Correlated Curve Trend (CCT) experiments were based on the concepts of O'Connor 



(1935 ). A more modern view is given by Bredenkamp (1984 ). In the CCT experiments, the 



desired stand densityn is achieved by thinning in advance of competition. An analysis of 
one of these experiments, based on the growth of Eucalyptus grandis (Hill) Maiden, is given 



by Bredenkamp and Burkhart (199C). In their paper they evaluate the use of various ways 
of quantifying the degree of crowding within a plot. 

Data from a CCT experiment established at the Border forest plantation in what is 
now known as Kwa-Zulu Natal, South Africa were used. In November 1936, test plots of 
Pinus roxburghii Sargent, a pine native to the Himalayas, were planted at an espacement 



of 1.80 X 1.80 m ( [Falkenhagen 19971 ). The area of each plot was 0.08 ha (800 m^). A 29 m 
wide buffer zone of trees was planted around each plot. The geographical details of the plots 
are given in Table ^.1\ 

Twenty measurements of the diameter at breast height^] (DBH) of each tree were taken, 
see Appendix H Measurements were usually taken during the summer months: October to 
March. At age 14, two measurements were taken, one in February and one in December. For 
this study these were averaged to give one measurement for that year. Height measurements 
were also made, but only diameter measurements will be examined in this report. 

In Figure Ej all the tree measurements for plot 2 are displayed. The mean of the 
measurements is also plotted. Each mean point is joined to its neighbouring mean points 
by piece-wise straight line segments. 

^Thinning is the artificial removal of trees by the forester. 

■^Tree stems per unit area. 

^The diameter at the breast height of the forester. 



Latitude (S) 


30°33' 


Longitude (E) 


29°45' 


Altitude 


1067 m 


Average annual rainfall 


945 mm 


Length of dry season^ 


3 months 


Mean annual temperature 


16.1°C 



■^Number of months with rainfall less than 30 mm. 



Table 2.1: Geographical information of Border forest plantation in Kwa-Zulu Natal, South 
Africa (After [Falkenhagen (1997| ).) 
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Figure 2.1: Individual growth curves for trees in plot 2. 



The average DBH vs. time curves will be modeled for each plot. Figures 2.2 and 2.3 
show the mean of the measured DBHs plotted against age for each plot used in this study. 
The number of trees in a plot at each measurement is also displayed as a vertical line at 
the age of each measurement. As can be seen in some of the plots, a particularly drastic 
thinning can cause a discontinuity in the growth curve. This is particularly noticeable in 
plot 15. 

In plot 1 the last few years of measurements actually show an increase in the rate of 



growth. This may be due to the natural decrease in trees within plot 1. Bredenkamp and 



Gregoire (1988) note a similar kind of behaviour in the Eucalyptus data which they studied. 



2.3 Growth Curve Modeling 

For modeling growth curve data, a large variety of models have been proposed. A review is 
given in Chapter 7 of 3eber and Wild (1989| ) . Approaches to growth curve modeling can be 



put into several different categories. 



2.3.1 Linear Approach 



This usually entails fitting polynomials to growth curves (Kshirsagar 1976). Some of the 



polynomial coefhcients can be assumed to be common for several different growth curves. 



Polynomial interpolation has been criticized as being biologically unreasonable ( Seber and 



Wild 1989). It is not commonly used in the forestry growth curve literature. 



2.3.2 Autocorrelated Errors 

Data which are collected from the same subjects at different times is known as longitudinal 
or time series data. Such data are often considered to have autocorrelated errors, (Eeber 



and Wild 1989). For example in modeling the growth in the weight of an animal, there 
might be a series of negative errors over a period of time that the animal was sick. For 
trees, correlated errors might be due to long periods of abnormal climatic conditions such 
as drought. 

Presumably, the autocorrelation in errors is going to depend on the frequency of mea- 
surements. If the measurements are far enough apart there is unlikely to be any autocorre- 
lation in errors. Also, stochastic analysis is easiest when there are equally spaced measure- 



ments, but there are ways of overcoming this restriction (McDill and Amateis 1991). 

In stochastic analysis, a difference function of the data is generally modeled by some 
differential form. In the next section, ways of modeling the data directly will be looked at. 
The basic differential forms the models are based on will also be discussed. 

2.3.3 Nonlinear Models 

Many nonlinear models have been proposed for growth curves. Some of them are based on 
biological principles. However, these biological motivations are not generally accepted as 
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Figure 2.2: Forest tree mean density at breast height growth samples for plots 1 to 10. The 
x's represent the mean DBH measurements and the vertical lines, the number of trees. 
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Figure 2.3: Forest tree mean density at breast height growth samples for plots 11 to 18. The 
x's represent the mean DBH measurements and the vertical lines, the number of trees. 



being compelling. Other models are purely empirical. The parameters of nonlinear growth 
curves can often be interpretable in terms of physical growth. 

Exponential and Monomolecular Growth Curves 

In the simplest of organisms, growth takes place by the binary splitting of cells. From which 
it follows that the rate of growth will be proportional to the current size of the organism, /: 

-r = nf, where k > 
dt 

which leads to the exponential growth curve: 

fit) = e'^(*-^) 

where 7 is a constant. Many growth models are exponential for small time, t. However 
an exponential growth curve model leads to unlimited growth, whereas growth is known to 
stabilize: 

a ~ lim f{t), 

t — *oo 

which implies 

— !■ U as r ^ 00. 

dt 

A simple way of achieving this is by assuming that the growth rate is proportional to the 
remaining size: 

rff 

— = K(a — /) with K > 0, 
dt 

which has the general solution: 

(2.1) f{t)^a~{a-f])e--'. 

If the function is used to describe monotonically increasing growth then 

a> (3>0. 
In this parameterization, a will be the final size, /3 = /(O) is the initial size and k governs 



the rate of growth. Equation 2A is generally known as the monomolecular growth model. 

Often in growth data, the growth first accelerates and then decelerates to a plateau. 
This gives rise to the "sigmoid" shaped growth curves. The point of inflection is the time 
when the growth rate is greatest. 

The logistic model has the following differential equation: 

df K 

-— — —f(a — /), where k > and < / < a. 

dt a 

Here the / in the numerator represents the tendency for the tree to grow indefinitely and 
the a — f represents the limiting component of the growth. [Zeide (1993 ) shows how most 
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growth equations can be broken up into an expansion and decline component. The general 
solution of the logistic model is 

a 



(2.2) 



/W- 



-oo < t < c». 



1 + e-'^(*-'r) ' 
The point of inflection occurs at t = 7 and the growth curve is symmetrical about this point. 



The Richards Model 

This restriction of the growth curve being symmetric about the point of inflection is not 
present in the Richards model ( Richards 1959| ). The differential equation for this model is 
given by: 

5-1 

- 1 



dt 



1-5 



I 



8+\, 



which leads to the model given by: 

This equation has enjoyed extensive use in forest tree growth modeling. However, it has 



also been subject to much criticism. Ratowsky (1983) shows that its parameter estimates 
are unstable. 

The Richards equation has an upper asymptote. Bredenkamp and Gregoire (198q ) 
note that tree diameter growth can start to increase again after reaching what appears to 
be an upper asymptote, due to tree mortality. 



The Schnute Model 

An equation which is more stable than the Richards model and also allows the possibility 



of non-asymptotic growth was introduced by Schnute (1981 ). Unlike most other growth 
models, the Schnute model is based on the acceleration of growth: 



(2.3) 



d'f f, ^k2\df 



The integrated version of the Schnute model is given by: 

-a(t-ri) nl/* 



(2.4) 
where 



fit) 



/r+(/2-/r) 



1 



1 — e-a(T2— Ti) 



t = age of interest 

Ti ~ youngest measured age 

T2 = oldest measured age 

/i = /(ri) 

/2 - /(T2) 

a = constant acceleration in growth rate 
b = incremental acceleration in growth rate. 



11 



In fitting the Schnutc model the initial guesses for /i and /2 are the initial and final measured 
size values. 

2.3.4 Estimation Methods 

Usually each growth measurement is assumed to have been drawn from a Normal distribu- 
tion: 

(2.5) y,\t,,9r^Noi-malifiU,9),a^). 

If data, {yi,ti} = {(yi, ti); (y2, ^2); ■ • ■ ; (2/Ar,iAf)}, are available then the maximum likelihood 
estimate (mle) of the parameters is given by: 

^mic = ma^p{{yi}\{ti},9) 

9 

N 

= maxexp - ^{Vi - f{U,d)f 
1=1 

N 

(2.6) =lmnY,{y^-f{t^,e)f 



i=l 



i.e. in this case the maximum likelihood parameter estimates arc found by minimizing the 
mean square error (MSE) of the predicted diameters by iterative numerical methods. 

Given 9, the errors, ei = yi — f{ti,9), are assumed to be independent. Thus, one 
diagnostic of the model fit is given by plotting the residuals and seeing if there is any 
unusual pattern of runs of positive and negative residuals, ( [Draper and Smith 198l| ). 

2.4 Modeling of Pinus Roxburghii Data 



i^'alkenhagen (1997) has studied nine different growth models for the diameter growth of the 
Pinus Roxburghii data discussed in Section E^. He found that the Schnute model, equation 



(2.4) had the least problems with convergence and gave the overall minimum mean square 



error. 



Zeide (1993) noted that the differential form of Schnute's model (equation (2.3)) does 



not provide a good fit for a Norway spruce tree data set. However, as the objective is to fit 



the integrated form of Schnute's model (2.4), this is not strictly relevant 



A good fit is obtained when all the data in each plot's data set is used ( Falkenhagen 



1997 ). Thus the Schnute model adequately interpolates the data. ( Falkenhagen 1997 ) also 
found the errors in the interpolation showed little autocorrelation (i.e. there were no unusu- 
ally long runs of positive and negative residuals), thus indicating that stochastic analysis 
(Section ^.3.2| ) may be unnecessary. A more difficult problem would be to see how well the 
model extrapolates the data. Thinning continued in some cases until the age of 24. It is 
pointless trying to extrapolate the tree growth while thinning is still taking place as then 
the model would have to predict the occurrence of any future thinning. As only age will 
be used as an explanatory variable, this would not be possible. Thus it will be required 
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that the model extrapolate from age 30 years onwards. By which time the tree has resumed 
normal growth. 

As will be seen in Chapter o, the Schnute model provides poor extrapolations for some 
plots and thus other modeling techniques need to be investigated. 
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Chapter 3 

Bayesian Statistics 



Probability, then, can be thought of as the mathematical language of uncertainty. 
R. L Winkler 

In this report Bayesian Statistics will be used in Neural Network Modeling, see Chapter 
y, and in comparing the performance of different regression techniques, see Chapter |5|. In 
this chapter all the Bayesian theory which is relevant to later work will be reviewed. Many 
aspects of Bayesian Statistics which will not be relevant to the rest of this report will not 
be discussed. Some of the results of examples in this chapter will be relevant to later 
developments. 



There are many text books on Bayesian Statistics, for example those written by Bergei 
(1985|), [Press (1989| ), |Box and Tiao (1992| ), [Bernardo and Smith (1994|), [Celman et al. (1995| ) 
and Jaynes (1996 ). Several, such as the book written by Jaynes (1996[ ), assume very little 



statistical training. The books by Box and Tiao (1992) and Jaynes (1996) are more oriented 
towards the natural sciences. 



3.1 Foundations 

Broadly speaking there are two schools of Statistics, Bayesian and Frequentist (also known 
as Classical or Orthodox). The Bayesian school is a minority, but has seen rapid growth in 
the last few decades. 

Frequentists generally only use probabilities to describe the proportion of times an 
event will occur in a given population. For example, if a rod is measured, a Frequentist will 
not consider its true length as a random variable. However, if there is a whole assembly 
line of rods then the length of the different rods within the assembly line can be assigned a 



random variable. A Classical view of the Bayesian / Frequentist debate is given by Papoulis 
|(T9§). 

Bayesians on the other hand use probability to express all types of uncertainty. A 
probability of zero corresponds to an impossible event and a probability of one to a certain 
event. Probabilities between zero and one express the degree of uncertainty. So in the 
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Bayesian framework it is possible to pose questions such as "What is the probabihty of a 
theory being true?" 

In a Bayesian sense, random variables are used to express uncertainty. Probabilities 
can be given for different proposed values of the variable. In Bayesian Statistics, model 
parameters are treated as random variables. 

3.2 The Rules of Probability 

A system for dealing with uncertainty that satisfies a certain number of reasonable desired 
properties, must be consistent with the following two rules ( laynes 1996| ): 



Product Rule : 

(3.1) P{AB\C) = P{A\BC)P{B\C) = P{B\AC)P{A\C). 

Sum Rule : 

(3.2) P{A\B) + P{A\B) ^ 1 

where A, B and C are propositions, e.g. 

A — A measurement of a quantity X will lie somewhere between xi and a;„. 

B = The samples {xi, a;2, ... , a;„} are measurements of X. 

C = A Gaussian probability distribution function should be used for X. 

The notation P{AB\CD) reads the probability of A and B being true given that C is 
true and D is false. A proposition is a statement that can be either true or false. Prior 
information will generally be denoted by an /. 

Since a continuous variable can take on an infinite number of values, its probability of 
being any particular number is infinitely small. Thus when dealing with random variables 
it is useful to work with a probability density function (pdf). The pdf for a variable X is 
defined as: 

(3.3) p{x\I) = ^P{X < x\I). 

Multivariate pdfs are defined in a similar way: 

p(x,y\I) = ——PiX < x,Y < y\I). 

The pdf, p{x\I), is just a function of x. Note that it could well have a different functional 
form to p{y\I). To distinguish the functional form, a subscript will be used. E.g. pxix'^\I), 
which is just the same function as p{x\I) except with all the x's replaced by x^'s. A random 
variable X has been distinguished from an instance of that variable, x. However, the same 
symbol will usually be used for the random variable and for an instance of that variable, 
but the meaning should be clear from the context. 
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It follows from the definition of pdfs, equation ( ^.3| ), that they must always be positive. 
The probability of a variable taking on a value contained in T) , which is a subset of the 
whole domain of the variable, V, is given by 



P{X eV')= f p{x) dx. 
JxeV 



From the sum rule, equation (p|2), it follows that 

P{X e V) + P{X V') = 1. 

From which it follows that 

(3.4) / p{x)dx = l. 

JxeV 

In the above, the probabilities have not been conditioned on any prior information, however 
all probabilities are based on some prior information and when it is not explicitly stated, it is 



assumed. Jaynes (1996 ) discussed the importance of bearing in mind the prior information 



a probability is based on. 

3.3 Bayes' Rule 

Using the product rule, equation ( |3.l[ ), Bayes' rule can easily be deduced: 

Many problems solved by Bayesian analysis take on the following form: 

1. Some data x = {xi, ... , a;„} are available. 

2. A pdf, p{x\9) is proposed, where 6 = {0i,02}, are a number of unknown model pa- 
rameters. This pdf is known as the likelihood. It contains all the information about 
how the parameters are related to the data. 

3. A prior pdf, p{9)^ which reflects the available prior information, is chosen. 
Writing Bayes rule in terms of the above pdfs, gives: 

(3.6) p{e\x) = -^g^- 

The pdf, p(Q\x') is known as the posterior distribution of 9. 

3.4 The Predictive Distribution 

The posterior pdf p(0|a;) must be normalized, i.e. 

p{9\x)d9 = 1. 
iVe 
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Using Bayes rule, equation (3.6), in the above equation, it follows that 

f Pimm ,, ^ , 

JeeVe Pi^) 
From which the marginal pdf, p{x) can be evaluated: 



(3.7) p{x) = / p{x\e)p{e)de. 

Note that p{x) depends on the form chosen for the likelihood and the prior. One should 
write p{x\I), where the / specifies the functional forms chosen for p{9\x,I) and p{9\I). So 
given different prior information (assumptions), /i and I2, the functional forms of p(a;|/i) 
and p{x\l2) can be very different. 

The pdf, p{x\I), is usually known as the predictive distribution, as it gives the prob- 
ability density function of any new measurement of data. When conditioned only on the 
prior information, p{x\I) is the prior pdf for the data. To get the pdf of some new data, 
Xn+i, given some old data {xi, . . . ,a;„}, the predictive distribution is as follows: 

(3.8) p{Xn+l\{xi,.. . ,Xn},I) = p{Xn+l\d,I)p{9\{xi,.. . ,Xn},I)d9. 

Je 



The pdf, p{x\I) is also referred to as the evidence, (see MacKay 1992a| ;p|). The reason for 



this is that it can be used to compare hypotheses. Say you have two hypotheses (theories) 
Hi and iJ2, and some data x. In order to compare the two hypotheses, in light of the data 
X — {xi, . . . , a;„}, the posterior odds ratio can be evaluated: 

PjHilx,!) _ p{x\Hi,I)P{Hi\I) 

P{H2\X,I) p{x\H2,I)P{H2\I)' 

So p{x\H,I) indicates how much the data contribute to the probability of H being true. 



i.e. what evidence it provides for H. An interesting example is given in Jcfferys and Berger 



(1992), where a fudged Newtonian theory and Einstein's Theory of General Relativity are 



compared in this manner. 



3.5 Eliminating Nuisance Parameters 



If one is interested in all the parameters, 9, then the posterior pdf of equation (3^) gives all 
the available information about 9 given the prior information and the data. For instance, 
the probability of the parameters being in a particular region T) is given by 



P{9 e V'\x,I) = / p{9\x,I)d9. 
JeeV 



This will not usually coincide with Frequentist confidence intervals. However, if only a 
portion of the parameters, 9i, are of interest and the rest of the parameters, 6*2, are nuisance 
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parameters, then the pdf for the parameters of interest can be obtained by 



p{ei\x,i)^ I piei,e2\x,i)de2 



(3.9) = / p{ei\e2,x,i)p{e2\x,i)de2. 

26^82 



This relationship follows from the general form of equation (3.7) 



3.6 Prior Probability Density Functions 

In this section, methods of choosing the prior p(9) are discussed. There are many ways of 
choosing the prior, and Berger (1985| ) has given a comprehensive review. Only those that 



will be relevant to the problems that will be solved in this research report will be looked at. 

3.6.1 Non-Informative Priors 

Often in scientific work, it is considered desirable not to include any information in the prior 
pdf, p{9). The interest is in what the data imply about p{9\x). Also there may simply be 
no useful prior information about the values of the parameters. 
From Bayes' rule, 

p(P\x) cc p{6)p{x\0). 

The posterior is only effected by the data via the likelihood function, p{x\6). So if one 
does not in any way want to distort the effect of the likelihood, the prior should be as flat 
as possible in the areas where the likelihood has an appreciable value and not have any 
relatively large fluctuations outside that area. For instance, if the prior had a peak in the 
tails of the likelihood, that could lead to an appreciable peak in the posterior. Also, if 
the prior was varying rapidly across the area where the likelihood had most of its mass, 
this would distort the shape of the likelihood. So, qualitatively speaking, non-informative 
priors should be as broad and featureless as possible. For a more detailed discussion on the 
qualities a prior should have see Berger (1985| ). 



In cases where little or no prior information will be used, many Orthodox Statisticians 
argue that Bayesian methods are not appropriate. In the Maximum Likelihood method, 
parameters are chosen which maximize the probability of the data, e.g. 

(3.10) ^mle ~ maxp(a;|0). 

6 

However, this is the same as making a m,aximum, a posterior (MAP) estimate and choosing 
a uniform prior for the parameters. The uniform prior is given by 

(3.11) p{0)=a 



where a is some constant. This prior cannot be normalized as in equation (3.4), it is therefore 



said to be improper. A uniform prior can be thought of as the limit of a Gaussian prior as 
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the variance goes to infinity. Improper priors can still lead to proper posteriors: 

p{x\e) 



(3.12) 



p{e\x) 



Jp{x\e)de 



where the constant, a, cancels out in the denominator and numerator. 



In the Bayesian approach, equation ( 3.12 ), the whole pdf of 9 is obtained. To obtain 
a point estimate, as in the Maximum Likelihood case, equation (p3n), a loss function is 
needed, see Section 



The Problem with Uniform Priors. 

One difficulty in assuming that p{6) is uniform, is that if some one to one function of 9 is of 
interest, = f{9), then p{(j)) will not in general be uniform. In order to determine the pdf 



of a function of a parameter, the following formula can be used ( Papoulis 1990 ) 



(3.13) 



P</.(</') ^Pe[<j)) 



df-'W 



where / ^{<p) denotes the inverse function of f{(f)). For clarity, subscripts are being used to 
distinguish the different functional forms of p(-). For example, if (f) — 9^^ then 

So, if pg{9) ~ 1, then p,p{(p) = <p~'^- Thus, by assuming complete uncertainty about where 9 
is, it is assumed that 9~^ is more likely to be closer to zero than further away. If, instead, 
d~^ was initially the parameter of interest, then the reverse would hold. Thus, the priors 
assigned to all the different functions of 6 are completely determined by which function of 
9 the uniform prior is assigned to. This is undesirable, since the choice of which function of 
9 to assign the uniform prior to is fairly arbitrary. 

Jeffreys' Priors 

As soon as a prior is assigned to some function of 9, this automatically implies what priors 



are assigned to all other one to one functions of 9, via equation ( 3.13 ). To make this whole 
family of priors invariant to which function of 9 was initially selected, the Jeffreys' prior 
can be used: 



(3.14) 



Pi9) (X [J{9)]'/\ 



where Jl{9) is the Fisher Information for 9 (Bernardo and Smith 1994): 

dHog{p{x\9)) 



J{0) 



(3.15) 



-E 



d9^ 

dHog{p{x\9)) 
d9^ 



-p{x\9) dx 
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where E[-\-] is the conditional expectation. Once the Jeffreys' prior has been assigned to 6, 
then any one to one function of 6, such as 4> — f{9) is also assigned a Jeffreys' prior: 



[pjWf (X J{c^) 



= -E 



= -E 



<P\ogp{x\4>) 


--f 




d02 

(f\ogp{x\e^ 


-\<p)) 


de 

d4 


2 





d6i2 







by defn. equation (3.14) 
by equation ( 3.15| ) 



J{0) 



de 



by the chain rule 
by equation ( 3.15| ) , 



from which it follows that 



pj{4>)^pj{0) 



de 



where the J" subscript is used to indicate that the prior was formed by Jeffreys' procedure, 
equation (3.14). As can be seen from equation ( 3.13| ), this is how pdfs should transform 
when a function of a parameter is taken. Thus, by choosing Jeffreys' prior for one function 
of e, all other functions of 9 are assigned their own Jeffreys' priors. 

Box and Tiao (1992) and Bernardo and Smith (1994| ) give further justifications for 
assigning Jeffreys' priors. 

As an example, if the likelihood is a Gaussian, p{x\e) = Normal(/^, cr^), and a is 



assumed known, then using equation ( 3.14 ) it can be seen that the Jeffreys' prior for /i 
is p(/i) = 1. If instead the mean, /i, is assumed known then the Jeffreys' prior for a is 

Jeffreys' prior can also be extended to multi-parameter problems: 

1/2 



(3.16) 



Pie) = \j\ 



where \J\ is the determinant of the Fisher information matrix defined by: 

d'^p{e\x) 



(3.17) 



A 



de,de. 



■ dx, 



where 0^ is the ith parameter of the vector of parameters 9. 

There are cases when Jeffreys' priors do not give good results for multi-parameter 
models. For instance in the case of a Gaussian distribution, p{x\e) = Normal(/i, cr'^), where 
both /i and a are unknown, the Jeffreys' prior is p{^,a) = l/cr^. This leads to a posterior 



with undesirable properties, as shown on pg. 361 of the book by Bernardo and Smith (1994). 
One ad-hoc procedure that has been proposed for overcoming such problems, is to 
assume some of the parameters are a priori independent. So in the case of the prior for 
p(/j,, ct), for a normal likelihood: 



pi^l,a) 



^p{li\a)p{a) 
^ p(m)p(ct) 
= l/a 
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where the single parameter Jeffreys' priors, equation (3.14), are used for p{p) and p{a). This 



prior leads to a posterior with more acceptable properties. 



Another approach called reference priors has been developed, (Bernardo and Smith 



1994 ), which reduces to a Jeffreys' prior in the single continuous parameter case and does 
not have some of the problems associated with Jeffreys priors in the multi-parameter case. 
However, reference priors are beyond the scope of this report. 

3.6.2 Conjugate Priors 

In general, the posterior, p{9\x), and the evidence, p{x), are not easy to evaluate. Thus, it 
can be desirable to choose the prior, p{0), such that the necessary calculations will be made 
easier. Usually any prior knowledge that is available is of a vague form and so the form 
of the prior pdf is fairly arbitrary, provided its properties are consistent with the available 
prior knowledge. 

One suggestion is to find a prior pdf which when combined with the likelihood function 
leads to a posterior pdf with the same form as the prior pdf. These are called conjugate 
priors. They have the added advantage that they lead to a more interpretable posterior. 

Example 3.1 Consider the case where x is normally distributed with a known mean /i and 
a standard deviation of a, i.e. x ^ Normal(/i, cr^). Instead of working with the standard 
deviation, the precision r = 1/cr^ will be used. It doesn't matter what function of the 
parameter is used because one can always transform back to the function of interest using 
equation ( ^.13| ). If x consists of N measurements, then the likelihood is given by: 

p{x\n,a) ^p{xi,X2,. .. ,XN\^J-,a■) 

N 



(3.18) =llp{x,\t,,a) 



(2w2)-^/%xp Uy(x.-A*)V(2a2) 



where in equation ( |3.18 ) it is assumed that each measurement is independent of every other 



measurement, given that /i and a are known. Expressing the likelihood in terms of the 
precision, r, gives: 



p{x\ij,t)(xt ' exp I —ry ( 



(3.19) ex r^/2 gxp (-tNs^/2) , 
where 

(3.20) s^ = 1 y (x, - fif 

i 

is the sample variance. Note that any constants of proportionality in the likelihood are not 
necessary because the posterior will be normalized at a later stage anyway. If the functional 
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form of the pdf in equation ( 3.19| ) is looked at with r as the variable and everything else 



as constants, then it is a Gamma distribution in r. Assume the prior for r is a Gamma 
distribution: 

(3.21) pir\a,c.) =. i^^^^W^^i exp(-™/(2c.))/(o,oo)(r), 

where 

1, re (0, oo) 



ensures the probability is only non-zero for positive values of t. The parameters a and uj 
must satisfy a > and lu > 0. The mean of the Gamma distribution is given by 

(3.23) E{T\a,uj)=uj 
and the variance is given by 

(3.24) y ar {T\a, Lu) =2lj^/ a. 

The parameters a and lu can be chosen to correspond to the desired prior mean and variance. 



Using equation (3.21) for the prior, the posterior of r is given by: 



(3.25) p{T\x,fi,a,uj) ocr(^+"-2)/2exp(-iT(a/cj + iVs2) j 



The posterior is also a Gamma pdf. Thus equation ( 3.21 ) is a conjugate prior to the 



likelihood given in equation (^.19). One advantage of conjugate priors is that it is easy to 



interpret the role played by the prior in the posterior. As can be seen in equation (3.25), 
a plays the same role as N and l/to plays the same role as s^. Thus one could interpret 
the prior, p(T\a,uj), as being equivalent to a measurements which have a sample variance 
of l/w. From which it follows that a natural interpretation for 1/uj is the prior variance of 



In Section 4.7 of Chapter M a Gamma prior will be used for the precision of some model 



parameters and model noise. 

3.6.3 Empirical Bayes 

In empirical Bayes methods the data are used in estimating the prior. In Section 3.4 it was 
shown how the marginal distribution of the data, p{x\H,I), contributes to the probability 
of the hypothesis being true: 

p{H\xJ)^p{H\I)p{x\H,I). 

The hypothesis, H , consists of two components: the likelihood, £, and the prior, tt. The 
type II maximum likelihood (ML-II) prior is obtained by maximizing the likelihood of the 
prior: 

TT — maxp(a;|7r,£, /). 
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Usually the maximization is done over some restricted family of priors. Using ML-II priors 
violates Bayes' rule since the prior is no longer independent of the data. However they can 



be used as an approximation to a true Bayesian approach. Berger (1985 ) elaborates on this 
distinction. 

Example 3.2 Given the prior, p{9\t) — Normal(0, l/r), then a suitable value of r has to 
be chosen. Using the ML-II method, 

T — maxp(a;|T, /) 

r 

(3.26) ^ max I p{x\0,I)p{e\T, I) de. 

^ Je 

Parameters, like r, which determine the prior distribution arc known as hyperparameters. 

One of the problems with the ML-II method is that it docs not acknowledge any 
uncertainty there may be in choosing the hyperparameters of a prior distribution. In the 
next section it will be shown how this additional uncertainty can be included. 

3.6.4 Hierarchical Priors 

Given the functional form for a prior distribution, p{9\t,I) without known values for the 
hyperparameters, r, then a prior can be assigned to t which reflects any uncertainty about 
its value. 



Example 3.3 Using the same prior as in Example 3.2, but instead of assigning the ML-II 



value for r, a prior is assigned to r. Choosing a conjugate prior gives 

p{t\I) = Gamma(a,a;). 
Values now have to be assigned for a and w which reflect the prior uncertainty in r. 

The following derivation shows the relationship between the hierarchical and ML-II 
priors: 

p{e\x) = f p{e,T\x)dT 



p{0\t, x)p{t\x) dr 
(3.27) ^fp(^o\r,xfJrMx]rl^^_ 

Jr PW 

From which can be seen that the ML-II method is equivalent to making the following 
approximation: 

P{x\t) «C(5(T-T^^lg), 

where t -i^ is the maximum likelihood estimate of r and c = p{t ■i^)/p{x) is a scaling factor. 
Thus the ML-II method only uses the maximum of the likelihood, p{x\t), while hierarchical 
priors make use of the whole likelihood function. 
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It is always possible to integrate out the hierarchical prior to get a single level prior: 

p{e) = I p{9\t)p{t) dr. 



One advantage of using hierarchical priors is that they generally are equivalent to single 
level priors which have very flat tails. This means they are robust, i.e. the final posterior 
does not depend strongly on the precise form, e.g. the mean and variance, of the hyperprior 



(Berger 1985). So the final result should not be too sensitive to the values chosen for a and 



LU in Example |3.3| . 

3.6.5 Exchangeable Parameter Priors 

Another use for hierarchical priors is when one has several parameters which are exchange- 
able. By exchangeable it is meant that one has no prior knowledge for distinguishing or 
grouping one or more of the parameters from the others. Probabilistically, this can be rep- 
resented as the prior, p{9i, 62, ■ ■ ■ ^On), being invariant to permutations of the parameters 



(Gelman et al. 1995). The simplest form of an exchangeable distribution is to have each 
parameter, 9i, independently and identically drawn from a distribution which has hyperpa- 
rameters r: 

N 

(3.28) p{e\T)^l[p{e,\T). 

In general the hyperparameter r is not known and so it is necessary to integrate over the 
uncertainty in t: 

N 
P{0)^ / \{p{e^T)p{T)dT, 

where pir) is a hierarchical prior. After integrating out the hyperparameter r, the param- 
eters in p(&) will not in general be independent, i.e. 

p{d,\Bu... ,0,_i,0,+i,... ,0^^)^p(0,)■ 
This provides a good remedy to the problem of overfitting, which will be discussed in Section 



4J of Chapter \^ 

Example 3.4 Consider the weighted sum, 

y(t) = ^u,/,(t) + e 

i 

where e is the component of y(i) which is not explained by i, otherwise known as the noise. 
If the prior values of /i(t) are independent of i then the weights Ui can be modeled as 
exchangeable. Therefore the posterior pdf is given by: 

v{Q\{vA) ^ Jpmy,t},r)p{T\{y,t})dT. 

Since r will determine the dependence between the weights, and r will be determined by the 
data through p{T\{y,t}), it follows that the dependence between the weights is determined 
by the data. 
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3.7 Loss Functions 

Generally a Bayesian analysis will result in a posterior pdf, either for a parameter of interest 
or for a future data value. It may be desirable to just report a single guess for the parameter 
rather than the whole posterior pdf. In which case it is necessary to decide which value to 



report. There is an extensive Bayesian theory on how to make these decisions ( Berger 1985 ). 
In Bayesian decision theory a loss is associated with any decision. For instance, if one 
is trying to guess the value of some parameter 9, then the loss associated with a guess, 9* , is 
a function, l{9, 9*). There are many possible choices for Z(-, •) depending on the application. 
A common choice is the square error loss function: 

1(9,9*) = {9-9*f. 

Another common choice is to use the absolute error. In Bayesian decision theory the opti- 
mum decision is given by choosing the value, 9* , which minimizes the expected loss: 

9* ^imivE\U9,9*)\. 
For the square error loss function: 

9* =TmnE\{9-9*f] 

= min f {9-9*fp{9)d9 
from which it follows: 

±.j[9-9*fp{e)d9 = Q 
{9-9*)p{9)d9 = 
(3.29) 9* ^ f 9p{9) d9. 



Thus when a square error loss function is used, the optimum value to choose is the posterior 
mean of the parameter. 

The zero-one loss function, is zero if the guess is correct, 9* = 9, and one otherwise. 
Its minimum expected loss is evaluated by setting 9* = inauXep{9\x), i.e. the maximum a 
posteriori value. It follows that the maximum likelihood method is a special case of Bayesian 
estimation, with uniform priors and a zero-one loss function. It seems advantageous that the 
Bayesian approach makes use of the whole likelihood function when making point estimates, 
while the Frequentist approach only uses the maximum of the likelihood function. 

3.8 Bayesian Computation 

As has been shown, many Bayesian calculations involve solving integrals, for example: 
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1. Obtaining posterior pdfs can involve integrating out nuisance parameters (see Section 



3^, e.g. 



(3.30) p(t|x) = / p(T,7|a;)d7, 

where 7 could be one or more nuisance parameters. 
2. To obtain point estimates, moments of functions of a parameter need to be found 



(Section 3.7), e.g 



(3.31) f*i9)= fi9)pie\x)de 

Je 

where f*{6) is the best point estimate for a function f{9), using the square error loss 
function. 

It often happens that these integrals are not analytically tractable. In such cases, numerical 
approximations have to be resorted to. If the integral is of low dimension then numeri- 
cal quadrature techniques can be used. However, for high dimension integrals, numerical 
quadrature is too time consuming due to the curse of dimensionality , i.e. the computation 



time increases exponentially with dimension (Evans and Swartz 1995). In order to approx- 



imate high dimensional integrals, numerical methods which make use of the probabilistic 
structure of the integrals are employed. 

3.8.1 Monte Carlo Integration 



Equation (3.31) is the mean value of the function f{9) with respect to the pdf p(6'|a::). One 
way of approximating a mean value, is to take samples from the distribution p{9\x) and 
then work out the sample mean of the function of interest, i.e. 



I 



/(%(01x)d0«l^/(0W) 



where the 6**^'^ are drawn from the pdf ^(6*10;). Here the superscript is being used to denote 
the sample number, e.g. 9^'^' is the second sample. As the number of samples, N, increases 
the more accurate this approximation will be. 



To approximate the integral in equation ( 3.3C ), samples can be drawn from p(r, 7|x). 
Each of these samples will contain values for r and 7. To get samples for p(r|a;), one just 
discards the 7 values. Although, this procedure does not give an analytical expression for 
p{t\x), the samples will allow any quantities such as moments, quantiles, etc. oi p{t\x) to 
be approximated. 

In order to employ these approximation methods, it is necessary to be able to draw 
samples from pdfs such as p{9\x). For simple distributions, such as Normal and Gamma, 
there are standard routines for efficiently drawing samples. For example many computer 
programs have commands for generating univariate normal distributions (pelman et al 



1995|). However, for more complicated distributions, it is often necessary to resort to Markov 



Chain techniques. 
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3.8.2 Markov Chains 



Gilks et al. (1996) give a comprehensive treatment of Markov chains in Monte Carlo inte- 
gration. Markov chains are a sequence of values 0' where 

i.e. the pdf of any value in the Markov chain depends only on the previous value. 

Markov chains can be used to simulate the drawing of samples from a pdf. If a Markov 
chain can be constructed so as its values converge to samples from a pdf of interest, say 
p{0), then it can be used in Monte Carlo integration. By definition the values generated 
by a Markov Chain are not independent. This usually means more samples are needed to 



obtain the same accuracy as would be obtained with independent samples ( Neal 1996 ). 

Markov chains can take a number of iterations before they start to converge to the 
probability distribution of interest. Thus it is common practice to discard a certain number 
of the initial iterations. Deciding how many samples to take from a Markov chain is often a 



matter of practical expediency. A number of convergence criteria are given by Cowles and 



Carlin (1995|). Some common techniques for constructing Markov chains are now discussed. 



Gibbs Sampling 

Although it might not be possible to sample directly from a pdf, p{9), it might be possible 
to sample from a subset of 6 based on the rest of 9, i.e. draw 9i from p{di\0-i), where d-i is 
the set of parameters 9 without the subset 9i. If the parameters are split up into s subsets, 
then the Gibbs sampling algorithm proceeds as follows: 



|(i+i)) 



1. 


Draw6'J'+^^ irom p{9i\9^^{). 


2. 


Draw 61^'+^^ bom p{92\9^i\^_2^t 


3. 




4. 


Draw6ii'+^^ from p(6l,|6lL'^^^). 


5. 


Let 0(^+1) = {0f+^\ .. . ,0i:'+^^ 


6. 


Let i = 2 + 1. 


7. 


Goto 1. 



Then, the draws of 9i can be considered approximate draws from p{9). 

It may happen that it is not possible to draw from any of the conditional distributions, 
in which case Gibbs sampling cannot be used. 

The Metropolis Algorithm 

In the Metropolis algorithm there is a, proposal distribution, pt{9^''^^'> \9^'''>) . The pdfpt(6'*^*+^^|6'*^*)) 
does not necessarily have to have any relationship with p{9), but must satisfy 
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One example could be a Gaussian distribution whose mean is centered on 9^ . 
An iteration of the Metropolis algorithm proceeds as follows: 

1. Draw ^ from ptC^^^+^^lfi*^*'). 

2. Set 

n{i+i) _ I ^ with probability min(p(^)/p(6'), 1) 
I 6^^'^ otherwise. 

The values of 9^"^' will then be approximate samples from p{9). 

If the proposal distribution is centered on 6'^*^ then the Metropolis algorithm can be 
seen as proposing a new value 0(*+i) = 9^'"' + e, where e is a random vector in the space of 
9. Thus the values of 6*^'^ will follow a random walk. 

An analogy can be made by considering the p{9) to be a surface where the areas of 
high probability are low and those of low probability high on the surface. If the position of 
a ball on the surface represents 9 then the Metropolis algorithm can be seen as randomly 
shaking the surface. Sometimes the shakes will move the ball uphill but usually downhill 
towards the areas of high probability. The position of the ball at regular intervals can then 
represent the samples of 9. 

This random walk behaviour can make the Metropolis algorithm inefficient, especially 
if some of the parameters in 9 are correlated. 
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Chapter 4 



Artificial Neural Networks 



4.1 Introduction 



The name artificial neural networks (NNs) covers a broad range of computational methods. 
There is a whole branch of the subject, which tries to model real biological neural networks, 
which will not be discussed in this report. A common feature of artificial NNs is that they 
consist of many interconnected simple processing units. The basic philosophy behind many 
artificial NNs is to use an algorithm that mimics the methods of information processing of 
a biological NN . 

NNs are usually used in classification problems. Other applications include regression. 



such as time series modeling (Weigend et al. 1991), and control (Miller III et al. 1990) 



This report will be concerned only with feed forward multilayer perceptron artificial 
NNs which are, arguably, the most popular type of NN. 

There is a wide range of literature in the NN field. The influential historical texts 
include Minsky and Papart (196£; 1990| ) and Rumelhart et al. (1986 ). A good introductory 



exposition is given by [Haykin (1994 ). A more advanced treatment can be found in Kung 



(1993| ). Statistical perspectives on NNs can be found in pliplcy (1994| ) and phcng and 



Titterington (1994) 
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4.2 Multilayer Perceptron Structure 



The multilayer perceptron neural network consists of layers of neurons. Figure 4.1 shows 
a graphical representation. The first layer is known as the input layer and the last as the 
output layer. The other layers are known as hidden layers. The NN in Figure 4.1 has only 
one hidden layer. This is the most common choice. 

Neurons are connected from the left layers to the right layers. The number of neurons 
in the input and output layers are dictated by the function being modeled. The number 
of hidden layers and the number of neurons in each hidden layer is a choice made by the 
modeler. The general function for a one hidden layer perceptron NN is given by 

(Nh / p \ ^ 

y^ Vijcj)'^ ^ UjkXk + aj\ + ^ Wiixi + hi 
]=i \k=i / I 

The meaning of symbols in this equation are: 

X : A p dimensional input into the NN. 

The fcth dimension of x is denoted by Xk ■ 

The output of output neuron i. 

The activation function of the output neurons. 

The activation function of hidden layer neurons. 

The number of hidden layer neurons. 

The weight of the connection from hidden layer neuron j 

to output layer neuron i. 

The weight of the connection from input layer neuron k to 

hidden layer neuron j . 

The weight of the connection from input layer neuron I to 

output layer neuron i. 

The weight of the connection from the hidden layer 

bias neuron to hidden layer neuron k. 

The weight of the connection from the output layer 

bias neuron to output layer neuron i. 



hi-) 



V. 



(4.2) 



Ujk 
Wil 

ak 



In Figure 4.1 the bias neurons are represented as squares, they are like input neurons 



with a constant +1 input. The NN in Figure 4.1 does not have any input to output layer 
connections. Usually the hidden layer activation functions are chosen to be sigmoid logistic 
or equivelantly, in terms of function approximation abilities, hyperbolic tangent functions. 
The output layer activation functions are generally chosen to be the identity functions in 
the case of function approximation. The input to output weights are often fixed at zero, i.e. 
they are deleted. When the output is one dimensional, only one output neuron is required. 



This leads to a subset of the family of functions in equation (4.1) 



(4.3) 



Nh / P 

f{x) = ^ Vj tanh ^ Ujk.Xk 
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Figure 4.1: A feed forward multilayer perceptron artificial neural network (NN). It has two 
input layer neurons, one hidden layer with three neurons and an output layer with two neurons. 
The squares represent the bias neurons. See the accompanying text for an explanation of the 
notation. 



where unnecessary subscripts have been dropped. The tanh function is hnearly related to 
the sigmoid logistic function, so in terms of function approximation it does not matter which 



is used. Neal (1996 ) prefers the numerical properties of the tanh parameterization. For the 
rest of the report, the family represented by equation ( |4.3| ) will be used. Any extensions to 
multidimensional output NNs is usually straightforward. 



4.3 Nonlinear Regression 

The purpose of NNs in function approximation is to approximate some nonlinear mapping. 



g{x). The range of functions which can be approximated by NNs was examined by Cybenko 



(1989 ) . He showed that one hidden layer NNs can approximate any continuous multivariate 
function with support in the unit hypercube. That is with sufficiently many hidden layers, 
the error of approximation can be made arbitrarily small. Cybenko's result only provides 
an existence proof. It does not specify how to construct the network, nor how the error of 
approximation is related to the number of hidden layer units. 

When the function to be approximated has noise added to it, then the problem becomes 
one of nonlinear regression. For an additive noise model, the corrupted function, y{x), is 
given by 

y{x) = g{x) + e 

where g{x) is the original function and e is a random noise component. The NN needs to 
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approximate g{x) given y{x). 

As can be seen from equation (4.1), NNs are just a family of nonlinear equations. 
The architecture of the NN determines the precise form of the equation used. The main 
user adjustable choice is the number of hidden layer units, Nh. NNs are often referred to 



as 'nonparametric' as the values of the weights in equation (4.1) generally contain little 
interpretablc information. 

4.4 Multilayer Perceptron Training 

In order to determine the weights of the NN, training data are required. Denote N samples 
of the input, output variables by 

Z? = (x«,y«), z = l,...,7V 

were a superscript is being used to denote the data sample number so as not to get confused 
with the data dimension, which is denoted by a subscript in this chapter. Then the weights 

9 = (u, V, a, b) 

= {ujk,Vj,aj,b) i = l...Nh,k^l...p 

can be inferred. The weights are usually chosen by least squares: 

(4.4) r =min £:(£>, 61), 

where 9* are the optimum weights and £{D, 9) is the sum of squares error (SSE): 

N 

(4.5) £{D^6)=Y.^f{x(^,9)~y(^f. 

4=1 

This is equivalent to a maximum likelihood solution when the noise is assumed to be iden- 
tically and independently drawn from a zero mean Gaussian distribution. From a Bayesian 
perspective, it is a maximum a posteriori estimation with uniform priors for 9. 

Usually, standard gradient descent is used to perform the minimization in equation 



(4.4). The weights are randomly initialized and then updated 



(4.6) aF+^) = 0p' + 5 -^£{D, 



e=eu) 



where 5 is known as the learning rate and is chosen by the user. The jth iterations weights 
are denoted by ^(■'' and the ith component of the weights is denoted by 9i. The derivatives 



of the weights in equation (4^) can be recursively calculated from output to input by using 



the chain rule, a procedure known as hack propagation (Rumelhart et al. 1986). 

The weights are updated until some convergence criteria are met. The SSE, the rate 
of change in the SSE and the number of iterations can all be used. Other techniques for 



deciding when to stop training will be discussed in Section 4.6 of this chapter 
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Gradient based techniques, such as equation (4.6), can be prone to being trapped in 



local minima. Simulated annealing can be used to circumvent this problem ( Kirpatrick et al 



1983) 



Another technique which is used in order to try avoid local minima and to try to 



decrease training time is to add a momentum term to equation (4.(;) 



lU+i) 



gO) 



^j-^(^' 



e=e'-j) 



" |-^(^' 



where a is a constant set by the user. The consequences of adding momentum have been 



investigated by phansalkar and Sastry (1994| ). 

When there is a large number of redundant training samples, the pattern update method 
can speed up convergence: 



gO + l) 



30) 



^|-(/(."r 



«)-»">)' 



e=eu) 



As the iterations progress, the training samples are cycled through, with fc = 1, . . . ,N. See 



Haykin (1994) for further discussion. 



There have been many other techniques for improving the time it takes to train a 



NN ( Kung 1993 ). Second order methods such as conjugate gradients have been found to 
decrease convergence time (Ripley 1994). However, Saarinen et al. (1993| ) show that the 
Jacobian matrix of £ with respect to the NN weights is generally rank deficient making the 
NN training numerically ill-conditioned. This can be the cause of the long training times 
that are usually experienced. 
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4.5 The Bias/Variance Dilemma 

In NN modeling a choice of the number of hidden units has to be made. The more hidden 
units used, the better the NN wiU be able to fit the data. However, it has been found (see 



Haykin (1994 ) for an example) that the generalization ability of NN often starts to worsen 
once too many hidden layer neurons are added. This phenomenon is known as overfitting. 
It is sometimes described as fitting some of the noise as well as the signal. 

Geman et al. (1992) shows how any regression estimate can be broken up into a bias 
and variance component. Some of their results are summarized below. The NN function's 
output, f{x), will depend on the possibly multivariate input, x. The true output, y, shall 
be considered univariate. It will be assumed that the data samples (x, y) are drawn from a 
multivariate distribution, p{x,y). 

A reasonable measure of how well the NN predicts y, is the expected square error for 
a fixed x: 

E[{y-f{x)f\x]^ fiy-fix,D))My\x)dy 

where the integral is taken over all possible values of y. The expected square error can be 
decomposed as follows: 

E[{y - f{x)f\x] = E[{{y - E[y\ x]) + {E[y\ x] - f{x))f\x] 
= E[{y-E[y\x]f\x] + {E[y\x]-f{x)f 
,^ ^. + 2i?[(y - E[y\ x])\ x] ■ iE[y\ x] - f{x)) 

= E[{y-E[y\x]r\x] + {E[y\x]-f{x)r 

+ 2iE[y\x]-E[y\x])-iE[y\x]~fix)) 

= E[iy^E[y\x]f\x] + iE[y\x]-f{x))^ 

The first term of the sum, E[{y — E[y\ x])'^\ x], is the variance of y for a particular x and 
is unrelated to the NN prediction. Thus, the second term is the appropriate measure to 
evaluate the NN performance. The form of the NN function will depend on the training 
data, D. This shall be indicated specifically by writing the function as f{x,D). Only for 
some training data sets will the NN give a good approximation of y. To get a training data 
set independent evaluation of how good f{x, D) is as an approximater, the expectation over 
all possible training data sets (of a particular size) of the squared error can be examined: 

ED[{I{x,D)-E[y\x]f]. 



34 



The bias / variance decomposition due to Gcman et al. (1992) is as follows: 

ED[if{x,D)-E[y\x]f] 

= ED[{{f{x,D)-ED[f{x,D)]) + [Ed[I[x,D)] - E[v\x])f] 

- Ed[{I{x,D) - ED[f{x,D)]f]+ Ed[{Ed[I{x,D)]- E[y\x]f] 
+2ED[{f{x, D) - Eoifix, DmEoifix, D)] - E[y\ x])] 

- EdUHx, D) - Ed[I{x, D)]f] + {Ed[![x, D)] - E[y\ xf 
+2ED[f{x, D) - ED[f{x, D)]] ■ {Eoifix, D)] - E[y\ x]) 

= iED[f{x,D)]-E[y\x])^ bias 

+ED[{f{x,D) ~ ED[f{x,D)]f] variance. 

If the expectation of the NN prediction is different from the expectation of y given x then 
it is said to be biased. An unbiased function may still have a large mean squared error by 
having a high variance, i.e. by being very sensitive to the training data. 

By increasing the number of hidden layer neurons, the bias is generally decreased 
while the variance is increased. Therefore, choosing the number of hidden layers is a trade- 
off between increasing the variance and decreasing the bias. The variance can be decreased 
by introducing more training samples. Thus, the more training samples available, the more 
hidden layer neurons can be introduced to decrease the bias. 

4.6 Methods of Avoiding Overfitting 

A method of determining the optimum number of hidden layers to use is to partition the 
data into a test and training set. NNs with different numbers of hidden units are trained on 
the training data. The error that each NN makes on the test set is evaluated. The number 
of hidden units that gave the smallest error is then taken as the correct choice. The NN 
with that number of hidden layer units can then be retrained on the whole data set. The 
drawback of this technique is that it makes inefficient use of the available data. Also, the 
user has to decide which data to set aside as a training set and which data to use as a test 
set. 

Automatic pruning techniques have also been developed for NNs ( Le Gun et al. 199C| , 



Hassibi and Stork 1993 ). Initially a large number of hidden units are chosen and then as 
training progresses an attempt is made to determine which hidden units are redundant and 
remove them. 

An alternative to finding the optimum number of hidden units is to choose a large 
number of hidden units and then regularize the solution. Two widely used techniques for 
doing this are weight decay and early stopping. 

In early stopping, a portion of the data is set aside as a test set. As training progresses 
the error on the test set is monitored. When the test set error reaches a minimum, train- 
ing is stopped. Although early stopping is quite commonly used, it has little theoretical 
justification. 

Weight decay on the other hand is a utilization of the general statistical procedure of 



regularization. Instead of minimizing the sum of squared errors (SSE), as in equation (4.5), 
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a regularization term is also included: 

N Ng 



(4.8) SiD, 9) = 5](/(x«, 0) - y«)2 + a^ 
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where the larger a is the more 'smoothing' is performed. Other regularization functions 
such as the first or second derivatives of the NN with respect to the weights can also be 
used. The optimum value of a can also be chosen by setting aside a test set. 

4.7 Bayesian Artificial Neural Networks 

In the previous sections a maximum likelihood solution to the learning problem was dis- 



cussed. Several Bayesian approaches have also been suggested ( Buntine and Weigend 1991 



MacKay 1992a, Sarle 1995|, Neal 1996). A good introductory overview is given by Bishop 



(1995) 



As discussed in Chapter 0, when predicting a new output, y^^'^^'^ given the N input 
output training set pairs D and input x^^'^ \ a probability distribution can be obtained by 
integrating over the model parameters: 



Using Bayes rule, the posterior of the parameters can be expressed in terms of its likelihood 
and prior: 

piy^''+'^\D,x<^^+^^) oc fpiy'^^+^'>\e,x'^^+^'>)piD\9)p{0)d0. 



To obtain a point estimate for j/*^^+^\ a loss function has to be assigned. As discussed 
in Section 3.7 of Chapter y, assigning a mean square error loss function leads to a point 



estimate given by the mean of the NN output: 

(4.9) y("+i) = / f{x'^^+'\e)p{e\ D) de. 



Similarly, an absolute error loss function leads to the point estimate being given by the 
median of p(y(^+^'| _D, a;'-^"'"^''). Credibility intervals can be obtained from percentiles of 

p{yiN+l)\D^^{N+l)Y 

4.7.1 Neural Network Priors 



MacKay (1992a) suggested zero mean Gaussian based priors for the network weights: 

p{e) = p{u)p{v)p{a)p{h) 
p 
= p{v)p{a)p{b) \\_p{uk) 



fe=i 
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where the notation is defined in equation (4.2) and the lack of a subscript denotes a group 
of weights. For example: 

Uk = {Ulk, ■ ■ ■ ,UN^k}- 

Groups of the weights share a common precision: 

Ujk '^ Normal(0,r„j^) 

Vj ^ Normal(0, Ty) 

Qj ^ Normal(0, Ta) 

b '- Normal(0, n) 

where j varies from 1 to A'^^ and k varies from 1 to p. The symbols r„j. , Ti, , Ta, ti, denote 
hyperparameters. The prior for the NN weights conditioned on the hyperparameters is given 
by: 

P N^ 
p{e\Tu,Ty,Ta,n) = Y[Y[p{Ujk\TuJp{Vj\Ty)p{aj\Ta)p{b\n). 
k=lj=l 

This grouping of the weights is usually justified by the need to account for different scalings 
in the output and input data variables. However, the grouping of the variables with the 
same hyperparameters follows directly from the principle of exchangeability, see Section 



3.6.5, For example, the hidden to output weights, Vj, all play the same role in equation 
(4.3) and so are a priori exchangeable. 

Neal (1996) has analyzed the relationship between the prior distribution on the net- 
work weights with the prior distribution on the network output. Some of his results are 
summarized below. 



From equation ( O ) it can be seen that the contribution of hidden unit j to the network 



function has the following properties: 

(4.10) E[vjhj{x)] ^ E[vj]E[hj{x)] 

where hj(x) is the output of hidden layer neuron j: 



p 



hj{x) = tanh qj + N^ UjiXi 



The factorization of the expectation is possible since Vj and hj{x) are a priori independent. 
The prior expectation of Vj is zero by definition and so 

E[v,h,ix)]=0. 

The variance of the contribution of hidden unit j is given by 

E[iv,h,{x))^]^E[v^]E[h'^{x)]. 

Define 

Vix) ^ E[h'^{x)]. 
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The limits of V{x) are given by the tanh function: 

v^(x)e[o,i] 



As can be seen from equation 4.3, the output of the NN is equal to the sum of the contribu- 
tions of the output's of the hidden units and the output bias unit. As the number of hidden 
units, Nh, becomes larger, the Central Limit Theorem can be invoked to get 

fix) - Normal(0, NnalVix) + al). 

Thus, in order for the NN to have a stable variance as the number of hidden units increase, 
the hidden to output weights have to be scaled: 

al = Wy/Nn 

where w^ is the maximum variance the hidden to output layer weights are chosen to con- 
tribute. With this rescaling it follows that: 

fix) - Normal(0, w„F(a:) -I- cr^). 

The prior variance of the output function will therefore remain stable as the number of 
hidden layer neurons increases. 

Neal (1996) argues that this rescaling will counteract the tendency of the NN to over 
fit the data as the number of hidden units increases. From which it follows the only limit 
on the number of hidden layer units should be dictated by computational constraints. 

Williams (1995) uses a maxim,um entropy approach to argue that a Laplace, rather 
than a Gaussian prior, should be used for the network weights. He shows how this can be 
used to implement a Bayesian pruning algorithm. 

4.7.2 Computational Techniques 



With the hyperparameter priors, the point estimation procedure of equation (4^) becomes 



(4.11) y(«+i) = / fix^^+'^\0)pi9\D,j)pi-t)ded-f 

where 7 = {Tu,Ty,Ta,Ti,,Tn} with r„ representing the precision of the noise added to the 
data. Although the noise precision is not strictly a hyperparameter, it is grouped with the 
hyperparameters as it is treated in a similar fashion. The integral required for the solution 



of the posterior predictive solution in equation (4.11) is difficult to solve and several different 
approaches have been proposed. 

MELximum Posterior Density 

Sarle (1995) advocates obtaining point estimates for the predictive NN result, equation 



(4.11), by maximizing the posterior probability of the weights and the hyperparameters. 



The hyperparameters are given slightly informative conjugate priors, sarle (1995) gives 
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simulation results on which the maximum posterior approach outperforms maximum likeli- 
hood and early stopping methods. 

This approach has the advantage of being very computationally efficient. However for 



small data sets, evaluating the full posterior as in equation (i.ff ) will probably produce 



better results as the mode may not be representative of the whole distribution. 
Gaussian Approximations 



The modes of the NN posterior can be approximated by Gaussians (Buntinc and Weigcnd 



1991, MacKay 1992a, Thodberg 1993|). The modes of the posterior of the network weights 



are found by optimization. Each node is approximated by a Gaussian whose covariance 
matrix is chosen to match the second derivatives of the log posterior probability at the 
mode. The posterior predictive distribution of the NN is found by the weighted sum of the 
integrals of each of the modes. 

In MacKay's approach the network hypcrparameters are estimated by ML-II type 
methods (see Section 3.6.3| of Chapter 0). Buntine and Weigend (1991) advocate analyti- 



cally marginalizing the hypcrparameters instead. MacKay (1994) argues that this produces 
less accurate results than the ML-II method as the modes of the marginalized posterior dis- 
tribution of the network weights can be quite unrepresentative of the posterior distribution 
of the weights as a whole. 

4.7.3 Markov Chain Monte Carlo Integration 

Neal (1996) advocates the use of Markov chain Monte Carlo (MCMC) techniques (see Section 



3.8.2 of Chapter m to solve the integral in equation (4.11). This has the advantage of not 



requiring any approximations to the parametric form of the posterior. 

In this scheme samples need to be generated from the posterior for the weights, p{9\ D). 
To do this, samples can be generated from the posterior of the weights and hypcrparameters. 



p{9,^\ D). The integral in equation (4.11) can then be approximated by: 

Ng 

i—l 

where 9^"^^ is the ith sample of weights and Ng is the total number of samples. 

The posterior for the weights and hypcrparameters is given by multiplying the prior 
by the likelihood: 

N 

p{9,j\D)^p{j)p{9\j)Y[p{y^'^\x^'\9,j). 

Conjugate priors are used for the hypcrparameters t„^, Tf,, r„ and Ta. A conjugate 
prior can also be given to the noise precision, r„. All the hyperparameters are precisions of 
the Gaussian distribution. Thus, their conjugate priors are given by Gamma distributions. 
For example, the conjugate prior of Ty is given by: 

Pi'^v) = — ^, — ^s — r"''/^"^exp(-r„ai,/2w„) 
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where the mean, w^, and shape parameter, «„, can be chosen by the user so as to give a 
suitably noninformative hyperprior. Each of the other groups of parameters are assigned 
their own hyperprior mean and shape parameter. 

Gibbs Sampling Updating of Hyperparameters 



In the scheme suggested by Neal (1996), the hyperparameters are updated by Gibbs samphng 



(see Section 3.8.2 of Chapter O.) The hkehhood of a hyperparameter depends only on its 



corresponding weights: 



(4.12) piv,,... ,v^Jt,) = (27r)-'=/V,^/2exp Lr„^z;f/2J . 

Thus, for a given group of weights, wi, . . . , vn^, the pdf of the hyperparameter, t„, condi- 
tioned on the weights is: 



p{Ty\Vl,... ,VNh,Oiy,U}y) CX r^"''+^'')/2 



^ exp -r„ a^/ujy + X! ^H /^ 



From this expression it can be seen that the prior for t^ can be interpreted as specifying a^ 
imaginary parameter values, whose average squared magnitude is l/ujy. Vague priors for t^ 
can be specified using small values of a„ . 

As before the noise for each data point is assumed to be drawn from an identical 
independently distributed zero mean Gaussian distribution with precision r„. The likelihood 
of the noise is given by 

(4.13) p{y\x,e,Tn) = (27r)-^/V„^/2exp Ur„ ^(y(^) - /(x^^), 0))V2 j . 

Using a conjugate Gamma prior, 

P{rn) = ^"/f"/°^% °/^-^exp(~r„a/2c.) 
r(a/2) 

the posterior is given by 

p(T„|i?,0)ar^"+")/2-iexpUr„(a/c. + ^(y(=)-/(a;(=),0))V2)j. 

Using these conditional posteriors, the weight hyperparameters and noise precision can be 
updated using Gibbs sampling. 

Hybrid Monte Carlo Updating of Network Weights 

The priors for the weights conditional on the hyperparameters are given by equations of the 
same form as equation ([4.12). The likelihood due to the training cases is given by equation 
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(4.13). The resulting minus log posterior is: 

P Nh Nh Nh 

- \0g{pi9\ D, 7)) OC ^ r„, ^ ul +TaJ2a^^+^vY.'■ 
k—1 i—1 i—1 j—1 



,2 



(4.14) 



N 



TnY^iV^'^ - I{x^'\e)f 



The form of equation (4.14) is similar to the weight decay error function of equation ([l.q). 
However, here the objective is to average over the distribution of weight decay coefficients 
and network parameters, instead of simply maximizing the posterior. Also, the user does 
not have to specify exact weight decay coefficients, but can instead specify a broad prior 
distribution. There is no need for a hold out set to evaluate the weight decay coefficients. 



Equation (4.14) is not amenable to Gibbs sampling as it is infeasible to sample from 



the conditional network weight posterior. Thus a Markov chain Monte Carlo (MCMC) type 



approach is more appropriate, see Section 3.8.2. 



However, the standard MCMC method can be very slow to converge when there are 
correlations between the parameters. Any proposed jumps which don't have a similar cor- 
relation structure will be likely to lead to improbable parameter values. To avoid this type 
of behaviour a method which takes into account these correlations needs to be formulated. 



Neal (1993; 1996) proposed the hybrid Monte Carlo method which was first formulated 
by Duane et al. (1987 ) in a Quantum Chromo Dynamics context. The hybrid MC method 
consists of two steps. First a dynamical simulation is performed and then the final result 
is accepted with a certain probability. The dynamical simulation part involves associating 
each network weight, 9i, with a particle coordinate in a fictitious physical system. With each 
coordinate there is also an associated momentum parameter pi. For the fictitious physical 
system a potential energy is defined: 

£{e)^~log{p{e\D,-t)). 

The corresponding momentum components contribute to the kinetic energy of the system: 

Nb 2 

I— 1 

where the rui are the associated mass components. The Hamiltonian of the system is given 
by: 

H{e,p)^E{e) + K{p). 

The coordinates are made to evolve according to the equations of Hamiltonian dynamics: 



de, 

dr 


dH 

dpi 


Pt 


dpi 
dr 


dH 


d£ 

de, 
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where r is the fictitious time. The dynamical simulation method is mathematically anal- 
ogous to a ball on a surface whose height is defined by —log{p{9\ D,j)). The particle will 
always "roll" back towards the valleys (posterior modes). 

In order to simulate the Hamiltonian dynamics a discretization procedure is used: 

p{r + e/2) =p,{t) - - ^^" 



9,(T + e)=0,(T)+e 



p,{T + e)=p,{T + e/2) 



2 86, 

Mr + e/2) 



mi 

ed£{e{T + e)) 
' 2 09, 



where e is the step length. One simulated dynamics iteration consists of L such steps. 

In continuous Hamiltonian dynamics, the Hamiltonian does not depend on r. However 
in the discrete case, H will not be the same at the start and end of the iteration. The new 
state is accepted with a probability 

min[l,exp(-iJ(6l^,/) + i/(6'\p*))] 

where 6' and p' are values of the parameters and momentum at the end of a dynamical 
iteration and 9^ and p* are the values at the beginning of the dynamical iteration. 

After each dynamical iteration, the hyperparameters are updated by Gibbs sampling 
and the momentum is updated by drawing from a multivariate Gaussian distribution. 

The hybrid MC method avoids the random walk behaviour of ordinary MCMC sam- 
pling and allows NN model parameters to be estimated in a reasonable amount of time. 



Many other implementation details and variations are discussed by Neal (1996 ). Free soft- 
ware implementing the above techniques is available from the URL 
http://www.cs.utoronto.ca/~radford. 

4.8 Conclusion 

Multilayer perceptron feed forward artificial neural networks (NN) provide a flexible non- 
parametric approach to nonlinear regression. The problem of overfitting can be solved by 
regularization. The amount of regularization can be chosen by cross validation or Bayesian 
techniques. The Bayesian technique can be implemented using maximum posterior tech- 
niques, Gaussian approximations and ML-II techniques or by MCMC sampling. A hybrid 
MC approach is required for practical computation times. 

The Bayesian approach provides a natural way of incorporating regularization without 
the need for a hold out data set. It also provides a whole distribution instead of just 
a point estimate and so credibility intervals can easily be generated. In Chapter M the 
MCMC Bayesian NN implementation will be used to extrapolate the forest tree growth 
data discussed in Chapter ||. 
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Chapter 5 

Methods of Comparison 



In this chapter a statistical methodology for comparing two fitting methods is discussed. 
Specifically methods for comparing the Schnute and artificial neural network (NN) extrap- 
olations of the forest tree growth data discussed in Chapter g are examined. 

Frequently when two methods are compared in the NN literature only the difference 
in the fits is given, e.g. the difference in the mean square error (MSE) on a test data set. 
However, it is also beneficial to determine how significant the observed difference is, i.e. 
whether or not the observed difference is only due to random variation. 

5.1 Criteria for Comparison 

As discussed in previous chapters, the parameters of a regression method are determined 
by training data. To test how well the method will predict other data drawn from the same 
distribution as the training data, a test set of data is usually employed. This is because 
the method will usually have an optimistically biased performance on the training set. It 
is possible for a method to be tailored to work very well on a particular data set. However 
this is no guarantee that the method will generalize well to other data sets drawn from 
the same distribution. For example, when using a polynomial of the same degree as the 
number of data points, the performance on the training set will be perfect but it is unlikely 
to generalize well. This phenomenon is known as overfitting the data. 

Following Rasmusscn (199(^ ), the factors that effect the evaluated performance of a 



regression method are defined as follows: 

1. The set of test data selected. 

2. The set of training data selected. 

3. The stochastic aspects of the method, e.g. random weight initialization and stochastic 
training.]^ 



'^Rasmussen distinguished stochastic prediction (e.g. Monte Carlo estimation) from the stochastic training 
element. 
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Thus an appropriate loss function (see Section 3.7 of Chapter 0) for a method trained on a 
training data set of size n would be: 

(5.1) Gpin) ^ / L[Fr^^rt{'^n,x),t]p{x,t)p{'Dn)p{ri)p{rt)dxdtdVndridrt. 



The functional form of the method is denoted by F. The L is the loss function for predictions 
made using training set 2?„, when the input is x and the correct test result is t. The r^ and 
rt denote the random initialization and training method. This is just the usual method of 
integrating out nuisance parameters that was discussed in Section ^^ of Chapter ^. 

In practice, Gp (n) can be approximated by averaging the loss over a large number of 
different experiments with different training and test sets. The dependence of n could also 
be removed by summing over experiments with different numbers of training cases. How- 
ever, the conditions of the experiments should be drawn from the probability distribution 
p(x, t, D, n, Ti, rt). In the next section a method of determining the statistical significance 
of the differences in G for two different methods is discussed. 



5.2 Hierarchical Analysis of Variance 



Kasmussen (1996) proposed the hierarchical analysis of variance (ANOVA) method for em- 
pirically comparing two regression methods. This technique will be applied in comparing 
the MSE of the Schnute and Bayesian NN methods. 

There are 1=11 forest tree plots available. Within each plot the last J = 4 points will 
be left out as a test set. The functions are each trained individually on each plot and their 
MSE is evaluated for each of the J test points. Let yij be the difference in the Schnute and 



NN MSE for test case j in plot i. Following Example 7 given by Spiegelhalter et al. (1996), 
the difference in residuals are modeled by 

yij r^ Normal(/i.,,rwithin) 
H^ ~ Normal(6',Tbotwoon)- 

The within plot variance, cr^jthin' ^^ ^^^ inverse of r^ithin- The between plot variance, 
'''between' ^^ ^^^ invcrse of Tbotwcon- The true mean difference between the techniques is given 
by 6. The cr^etween measures the difference between the plots. It can be interpreted as 
the difference caused by the different training sets in each plot. The cr^ithin measures the 
difference caused by estimating the MSEs from a finite number of test samples. 

The main interest is to evaluate if one method is significantly better than the other. 
How this is done from both a Frequentist and Bayesian perspective is now described. 

5.2.1 Frequentist Estimation 

An unbiased estimator of 9 is given by 

ij 
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In Frcqucntist terms one method is said to be significantly better if the p-value, given by 
pi& > |y.. I 1^ = 0), is less than some threshold. That is, the probability of having the current 
or more extreme data, assuming that the methods are the same on average, is examined. 
The smaller the p- value the more significant the observed differences are thought to be. 
The t statistic is given by 



t = y. 
where 



(7(7^ &-: 



y^■ ^j^vv- 



-1/2 



\2 



J 

3 



It has a t distribution with / — 1 degrees of freedom 

p{t) ^ fl + 



e ^-^/^ 



/- 1, 

From which it follows that the required p-value is given by 



t 



p = i-y p(t')dt', 

where the integral can be evaluated numerically using the incomplete beta distribution. 
Unbiased estimators are available for the variances: 



''within ^ Ji J — \\ ^-^ 2-^^y^^ y' 



(J-1) 
J 



''"between ^ W _ 1 A^^y^- y-'' '^within) I J- 

i 

These can be used to evaluate the cause of the variation. 

5.2.2 Bayesian Estimation 

In the Bayesian case there are two distinct ways of determining whether the two methods 
are producing significantly different results. The ratio of the hypothesis that ^ = (there 
is no difference in the true MSEs) against 6* 7^ can be examined: 

p{e ^ Q\y.,) p{e ^ 0) 

p(MO|y.,) p{Q + Qf^''''^ 
where B{yij) is the Bayes factor and is given by 

Assuming equal a priori probabilities for the two hypotheses, inference would be made using 
the Bayes factor. The smaller the Bayes factor, the more probable that the two methods 
are different. 
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The second Bayesian method would be to evaluate p{0 > 0\yij). A Bayesian treatment 
of hierarchical ANOVA is given in Chapter 5 of Box and Tiao (1994 ). They use the Jeffreys' 
prior: 

^\"j '''between' '^withiiJ ~ P(.^JPv'''between' ''^within/ 

with p{0) assumed uniform and 

/ 2 2 \ —2 —2 

^V^bctween' within/ between within' 

From which it follows that the posterior is given by: 

PiS\{yij}) -- a2"^'Betaia,/(ai+a2)(P2,Pl), 

where Betai is the incomplete beta distribution and 

i j 

I{J-l) 

The inferences drawn from Bayesian hierarchical ANOVA can be sensitive to the chosen 
prior. This is because the likelihood does not contain much information about Cbetween (^^^ 



\2 



a2 


=^)^iy 


P2 


I 

~2' 



page 18 and 19 of Elasmusscn (1996)) 



Alternatively, hierarchical priors could be used (see Example 7 of Spiegclhalter et al 



1996): 



9 '^ Normal(/ie, Te) 

Twithin ^ Gamma(awithin, /^within) 
'Hjetween ^^ '^^mi^^v^betweeni Pbetween j 



where the hyperparameters are chosen so as to be uninformative. In Rasmussen (1996) 
this approach is suggested but not implemented. Gibbs sampling needs to be used to make 
inference about the parameters. 

In the next Chapter, the hierarchical ANOVA technique will be applied in comparing 
the results of the NN and Schnute extrapolations of the forest tree growth data. 
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Chapter 6 



Results and Discussion 



In this chapter I present the resuhs of fitting the artificial neural network (NN) and Schnute 
functions to the forest tree growth data presented in Chapter g. Graphical representations 
of the NN and Schnute fits for each of the forest tree plots are given. Credibility intervals 
are also plotted for the NN fits. An analysis of variance (ANOVA) is used to compare the 
two methods. 



6.1 Experimental Design 

The main purpose is to determine which technique, the Schnute or Bayesian NN, is best at 
extrapolating forest tree growth data. As discussed in the previous chapter, in order to do 
this the results when different training sets are used need to be compared. A number of test 
cases are needed or else the variation due to sampling error, cr^^i^i^, will be too high. 

Eleven different sets of forest tree data were available. Each plot is extrapolated sepa- 
rately, i.e. for each plot the methods will be trained on a subset of the data and extrapolated 
on the remaining subset. 

Choosing how much data to use for training and how much to use for testing is a 
difficult problem. If too little data are used for training then the predictions will have a 
high variance. This was the role of o-bctwocn i'^ ^^^ ^^.st chapter. If too little data are used for 
testing then the sampling error will be to high. Ideally, data for the whole function should 
be tested, but for practical reasons only a small sample of the function's points can usually 
be checked. 

So in deciding how much of the data to divide into test and training sets, the effects 
on cr^ithin ^"^^ '''bctwcGn havc to be considered. As discussed in Chapter |[ there is a limit to 
the number of points which can be chosen as test data. The first 15 data points are affected 
by the removal of trees. As the number of trees is not included as an explanatory variable, 
this extra variation could not be taken into account by the prediction method. Therefore 
it was decided to restrict the data points used, for testing, to the last four. For some of 
the plots where thinning ended earlier, more data could be used. However, the statistical 
analysis of the results is easier when the number of test points is the same for each plot. 
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Extrapolation can be useful to the forest manager in making predictions for future 
growth. 

In summary, the first 15 points (ages 5 to 28 years) are used as training points and the 
last 4 points (ages 30 to 37) as test points. 

6.2 Implementation Details 

In this section I discuss specifically how the Bayesian neural network (BNN) solution and 
the Schnute solution were obtained for the forest tree growth extrapolation problem. 

6.2.1 Bayesian Neural Networks 

The NN will need one input layer neuron for the age and one output layer neuron for the 
average density at breast height (DBH). This study was restricted to one hidden layer NNs, 
as they are generally adequate for simple nonlinear regression problems. 

As discussed in Chapter 0, the number of hidden layer neurons can be chosen as high 
as computationally feasible in the Bayesian scheme, without the danger of overfitting. Eight 
hidden layer neurons were chosen, as for the relatively uncomplicated curves I wish to fit 
they should be more than adequate. The hidden layer neurons were given tanh activation 
functions and the output neurons were assigned linear activation functions. 

Radford Neal's Bayesian NN package, bnn, (see Chapter Q), was used to fit the data. 
In choosing the prior parameters and Markov chain Monte Carlo (MCMC) scheme, the 
regression example distributed with the bnn package was used. The choices appear to be 
fairly generally applicable. Although I did choose to normalize the training data to zero 
mean and unit variance, so as to be more compatible with Neal's regression example. 

The noise variance was given an inverse Gamma prior with a mean precision corre- 
sponding to a standard deviation of 0.05 and a shape parameter of 0.5. This corresponded 
to a prior distribution which spans several orders of magnitude and so should easily en- 
compass the noise level of the data. The hidden layer weights were assigned a zero mean 
Gaussian prior. The variance of the Gaussian prior was given an inverse Gamma prior with 
a scaled (see Chapter ^ mean equivalent to a standard deviation of 0.05 and a shape pa- 
rameter of 0.5. The input to hidden layer weights and the hidden layer biases were given 
equivalent priors but unsealed. The output layer bias was given a zero mean Gaussian with 
a standard deviation fixed at 100. Giving the output layer bias a hyperparameter was un- 
necessary because hyperparameters are only needed when there is more than one parameter 
in a group. 

The MCMC was split into two phases. In both phases, for every iteration, the momen- 
tum variables are updated by the heatbath method and the noise level is updated by Gibbs 
sampling. The other details of each phase are as follows: 

1. In the initial phase the objective was to get a fairly reasonable starting point. The 
hyperparameters were all fixed at the moderate value of 0.5 and the network param- 
eters to 0. Then one iteration was performed with the following specifications. The 
network parameters were updated using hybrid Monte Carlo with a trajectory of 100 
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leapfrog steps long, a stepsize adjustment factor of 0.2 and a window size of 10 (Neal 
19961) . The hyperparameters remained fixed. 



2. In the final phase the actual MCMC sampling is performed. For each iteration, the 
hyperparameters were updated using Gibbs sampling. The network parameters were 
updated using a hybrid Monte Carlo method with a trajectory of 1000 leapfrog steps, 
a window of size 10, and a stepsize adjustment factor of 0.4. 

It was decided to do 1000 sampling iterations for each plot, while in the bnn regression 
example only 100 iterations were performed. 

On a 66 MHz Cx486DX2-S PC with the Linux 2.0 operating system, each plot took 
approximately 50 minutes to train. Of course shorter times can be achieved by using less 
iterations, but this will effect the accuracy. 

Convergence 



For most of the plots, the Markov chains seemed quite homogeneous. Figure 3.1 shows 
a plot of the hidden to output layer hyperparameter for a NN trained on forest tree plot 



2. Although there are formal convergence tests ( Cowles and Carlin 1995 ) visually the time 



series looks fairly homogeneous. However there were some plots where this was not the case. 



Figure 6.2 shows the same hyperparameter for plot 1. Here the last 300 points have a much 



higher mean and variance than the first 700. Examining 5000 iterations. Figure 6.3 shows 
that the Markov chain seems to oscillate between these two 'states'. Adjusting the step size 
adjustment factor so as to modify the rejection rate did not seem to make much difference 
to this problem. The posterior may be multi-modal and the Markov chain is oscillating 
between the two modes. Examining the samples for the training MSE for plot 1, Figure 



6.4, shows that the MSE seems to be the same for both modes and so averaging over the 
different modes may still be valid. Taking the first 1000 iterations gave good results. Taking 
only the first 700 made no significant difference. 

Point Predictions and Credibility Intervals 

Point predictions were made by taking the mean of all the Markov chain (MC) samples 
for each point. As discussed previously, this implies a MSE loss function. The standard 
deviation for each predicted point was also evaluated. If the predictions were assumed to 
follow a Gaussian distribution then the true value to be predicted would have approximately 
95% probability of being contained within two standard deviations of the mean. 

The 95% credibility intervals of the predicted data can be estimated directly from the 
MC samples, without any assumptions about the distribution. I evaluate these for each 
predicted point as well. These intervals are distinct from the credibility intervals of the 
predicted mean in that they take into account the noise variance. The 95% credibility 
intervals were obtained using the same scheme as was used for the median evaluation in 
bnn. First the noise standard deviation cr„ is estimated by taking the mean of the Monte 
Carlo samples of the noise standard deviation. Then to predict the 95% credibility interval 
of the output, f{t) for an input t, each Monte Carlo estimate of f{t) is used to generate 200 
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Figure 6.1: Markov chain of hidden layer to output layer hyperparameter for plot 2. 
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Figure 6.2: Markov chain of hidden layer to output layer hyperparameter for plot 1. 
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Figure 6.3: Markov chain of hidden layer to output layer hyperparameter for plot 1. 
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Figure 6.4: Markov chain of average training MSE for plot 1. 
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samples from a Gaussian with mean f{t) and standard deviation ct„. Then I evaluate the 
value for which 2.5% of these samples lie below, and the value for which 2.5% of the samples 
lie above. Ideally the 95% interval should be chosen so as to be as short as possible, see 



page 85 of Box and Tiao (1992 ). However this would have been much more computationally 
expensive. 

Evaluating the standard deviation and 95% credibility intervals required modifying the 
net-pred command in the bnn package. 

6.2.2 Maximum Likelihood Neural Network Solution 



The maximum likelihood or least squares optimization for the NN fit (see Section 4.4 of 



Chapter m was also tested. The Nevprop packageF^ (Goodman 1996|) was used. The inputs 



and outputs were normalized to unit variance and zero mean. Ordinary back propagation 



performed better than many of the available alternatives, such as Quickprop, (Goodman 



199q ). The NN generally converged to the training data quite quickly, but many further 



iterations were required to obtain good test results. The test error appeared to slowly 
oscillate, with the oscillations gradually diminishing in amplitude. It was only after several 
hours of iterations that comparable results to the Bayesian learning solution were achieved. 
Increasing the learning rate lead to wild oscillations. Adding a momentum term was found 
to be of little help. 

Using a tanh output activation function improved the rate of convergence but required a 
proportion of the training data to be set aside to test when the optimal fit had been achieved 
(Gordon 1994a|jb[). Leaving out a validation set did not help with the linear activation 



function NN. Comparisons between Bayesian and least squares fits are discussed in Sarle 
|(1995| ), [Neal (1996| ) and [MacKay (1992a| ). Generally, the Bayesian solution is found to 
produce superior results. 

Confidence intervals can be produced for the least squares fit using linearization and 



other methods (Brittain and Haines 1997). It would be interesting to compare these with 



the credibility intervals produced by the Bayesian method. 



6.2.3 Schnute Implementation 



For the Schnute function implementation, the fits obtained in Falkenhagen (1997 ) were used 



To improve the rate of convergence to the least squares solution in fitting the parameters 



of equation (2.4), Falkenhagen set ti to the initial training age (5 years) and T2 to the final 
training age (28 years). The initial values for the parameters yi and y2 were set at the 
corresponding initial and final training average DBHs. 



The secant method (Draper and Smith 1981) was used to fit the parameters j/i, y2, a, 
b, separately for each plot. 



^Available from ftp://unssun.scs.unr.edu/pub/goodman/nevpropdir/. 
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Plot 


BNN 


Schnute 


1 


61 


177 


2 


1 


43 


4 


21 


25 


6 


3 


10 


8 


3 


1346 


9 


101 


25 


11 


13 


523 


12 


19 


11 


15 


19 


8 


16 


6 


86 


18 


61 


38 


Mean 


28 


208 


Std. Dev. 


32 


406 



Table 6.1: Table of test MSEs in units of mm^ 



6.3 Extrapolation Results 



The BNN and Schnute techniques were applied to 11 plots. Plot 10 was not included in the 
evaluation as the Schnute function failed to converge in that case. The results are displayed 
in Figures 6.5 to 6.15. The x symbols are used for the training data for each plot and the + 



symbols are used for the testing data. The Bayesian network prediction and 95% credibility 
intervals are displayed as solid lines. The Schnute prediction is displayed as a dashed line. 
For some plots the Schnute function became complex over t £ [0,40]. It is not plotted 
over these regions. The Bayesian credibility intervals appear to behave in the expected 
manner in that as the prediction gets further away from the data points the credibility 
intervals widen. This behaviour was also noted in the non-MCMC implementation described 



by MacKay (1992a). The credibility intervals are not symmetric as would be the case if only 



the standard deviation of the prediction was used to evaluate the credibility intervals. 

The displayed intervals may appear to be too wide. In theory, each data point should 
have a 95% probability of being contained within them. However, if the point (0, 0) is taken 
as valid, then in 2 out of the 11 plots (plot 6 and 8) the origin is not contained between the 
95% credibility intervals. This suggests in some sense the credibility intervals are not overly 
wide. 



The MSE test losses for both techniques are shown in Table 6.1. The results indicate 



that the BNN method has a smaller MSE and is less variable. However, as the Schnute 
method produces such variable results it is difficult to be sure that on average its MSE isn't 
roughly equal to the BNN MSE, i.e. the different mean MSEs may only be due to sampling 
error. To test this, the ANOVA techniques discussed in Chapter ^ were used. 

The results of the Frequentist and Bayesian ANOVA analysis are shown in Table p^ . 
The 'ConjugatC]^' column is the ANOVA result using Bayesian analysis with a conjugate 
prior (see Section 5.2.2 of Chapter |^.) This analysis was performed using the Bugs software 
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Frequentist 


Conjugate^ 


Conjugatea 


Jeffreys 


e 


180 


181 


182 


261 


'^within 


404 


439 


439 


- 


'^between 


364 


373 


372 


- 


/between C^^)) 


45 


38 


38 


- 


p{%) 


18 


8 


7 


11 



Table 6.2: Table of Frequentist and Bayesian ANOVA analysis. Where not specified, units are 



packageg. The mean was given a normal distribution with mean and precision 10^ ^*'. 
The between and within variance prior precisions were both given gamma distributions, see 
equation ( ^j.21 ), with a = 2 x 10~^° and w = 1 

The 'Conjugate2' column was generated in exactly the same way as the 'ConjugatC]^' 
column except that a was 1000 times smaller. This makes the distributions roughly 1000 
times more vague. 



The 'Jeffreys' column is the Jeffreys prior Bayesian analysis (see Section 5.2.2 of Chap- 
ter |5[) The cr^etween ^'^'^ ^within were uot analyzed for the Jeffreys prior analysis. 

The '/between' '^^^ ^^ ^^® percentage of the variance explained by the between variance 
component. The 'p' row is the p-value for the Frequentist analysis, with the null hypothesis 
being that the two methods had equal expected MSEs. For the Bayesian analysis the 'p' 
row contains the probability of the Schnute model having a smaller MSE than the BNN 
model. 

The main concern of the Bayesian ANOVA approach is sensitivity of the results to the 
prior chosen for cr^et^j^ecn- However, as can be seen the results appear to be fairly insensitive 
to making the priors more uninformative in the conjugate case. 

The p-value obtained is the probability of obtaining the current data or a more extreme 
version of it, if the two techniques were really identical under a mean square test error 
loss function. The obtained value of 0.18 does not meet the usual 0.05 significance level. 
Hence, from a Frequentist perspective, the hypothesis that the two techniques have identical 
expected test MSEs cannot be refuted. A contributing factor to this inconclusive result is 
the high level of variation in the Schnute results. The BNN has an estimated standard 
deviation an order of magnitude smaller than the Schnute method. This factor may also 
prove important when choosing a method to use. As discussed by Rasmusscn (199^ ), very 
large data sets are usually needed to obtain significant result from a Frequentist ANOVA. 

The estimation of the mean test MSE difference for the conjugate priors shows good 
agreement with the Frequentist estimate. Thus the conjugate prior analysis is taken as 
the Bayesian result. The discrepancy in the noninformative analysis confirms the negative 



opinion of the approach discussed in Section 5.2.2 . 

The Bayes result indicates that there is a 93% probability of the BNN approach having 
a smaller test MSE than the Schnute approach. It is difficult to directly compare the 

^Available from http://www.mrc-bsu.cam.ac.uk/bugs/mainpage.html 
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Bayesian and Frequentist results as they use very different criteria. It is my belief that the 
Bayesian approach provides a more direct criterion for comparing the hypotheses. To make 
a decision using Bayesian decision theory, the mean loss associated with each method is 
compared. The method with the smallest mean loss is the one which will be chosen. If any 
losses associated with computation time and storage are ignored, then the mean losses can 
be treated as the corresponding mean MSB's. Thus from a Bayesian perspective, the BNN 
procedure would be chosen. 



The Bayesian and Frequentist conclusions do not have to necessarily coincide. Berger 



(1985) gives several examples where they differ considerably. Thus the Bayes result can 
strongly indicate that the methods are different, while the Frequentist result may be incon- 
clusive. I believe this is true in the current case. 

Whether or not the techniques are significantly different is not the only relevant crite- 
ria. A practically insignificant difference can still be deemed statistically significant if there 
are enough samples. To determine whether the observed MSE difference of 180 mm^ is 
practically significant, it can be converted it into a basal area error (BAE) difference. Mul- 
tiplying by 7r/4 gives a BAE difference of 1.4 cm^. The average tree height for the study was 



about 29 meters for the testing period ( Falkenhagen 1997 ). The average stems per hectare 



was approximately 600. Treating the volume as basal area times height, an approximate 
improvement of 1.7 m"^ per hectare (100 m^) is achieved. A plantation can be thousands of 
hectares in size. 
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Figure 6.5: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 1. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.6: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 2. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.7: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 4. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.8: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 6. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.9: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 8. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.10: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 9. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.11: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 11. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 



62 



0.5 



0.4 



0.3 



03 



E 0.2 

Q 

0.1 



-0.1 



-0.2 




1 I r 

15 20 25 

age / years 



Figure 6.12: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 12. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.13: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 15. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.14: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 16. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Figure 6.15: ADBH BNN (solid lines) prediction with 95% credibility intervals for plot 18. The 
Schnute prediction is displayed as a dotted line. The x are the training data and the + are the 
testing data. 
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Chapter 7 



Conclusion 



In this research report I have apphed Bayesian artificial neural networks (BNNs) to the 
problem of extrapolating mean forest tree growth curves. I have also used an analysis of 
variance (ANOVA) to evaluate the BNN performance. 

In Chapter |l|, I showed how the functional form of a multilayer perceptron artificial 
neural network (NN) could be interpreted as a sum of logistic functions. As the logistic 
function is a common model for growth curves, the NN approach has some justification in 
terms of measuring the average growth of a plot of trees. 

Reviews of forest tree growth curve modeling and Bayesian statistics were given in 
Chapter g and g|. The emphasis was on providing the relevant context and background for 
the theoretical and empirical developments in the rest of the report. 

In Chapter y, a survey of the relevant artificial neural network theory was given. Par- 
ticular emphasis was placed on the Bayesian Markov Chain Monte Carlo (MCMC) approach. 
It was shown that an alternative way of justifying the prior assigned to the NN weights was 
to treat the weights as a priori exchangeable. Each group of weights that was assigned a 
common hyper-prior, played a priori equivalent roles in the NN functional form. Previously, 
the groupings in the NN prior were justified in terms of different input and output scalings. 

In Chapter g, analysis of variance (ANOVA) based techniques for comparing two re- 
gression methods were discussed. 

The main results of the research report were given in Chapter O. Eleven forest tree 
average density at breast height (DBH) growth curves were modeled. From each forest tree 
plot, nineteen age DBH pairs were available. Fifteen were used as training samples and four 
as test samples. The Schnute and BNN techniques were applied separately to each plot. 

The BNN technique had an average mean squared error (MSE) across the plots of 28 
mm^ with a standard deviation of 32 mm^ . The Schnute technique had an average MSE of 
208 mm^ with a standard deviation of 406 mm'^ . 

The Frequentist, Jeffreys and Conjugate prior ANOVA methods were applied. I do 
not know of any other instances of the Conjugate prior ANOVA being used in this context. 
The Frequentist and Conjugate ANOVA both yielded very similar parameter estimates. 
The Conjugate ANOVA did not appear to be sensitive to changes in the prior. The Jeffreys 
ANOVA parameter estimates were different, which probably lends support to the misgivings 
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other workers have expressed about this technique. 

The Frequentist ANOVA gave a p-value of 18%. This does not meet the usual 5% 
requirement. Thus, from a Frequentist perspective, the null hypothesis, that the two tech- 
niques will on average produce equal MSEs, cannot be rejected. This may be due to the 
small sample size and high variance of the Schnute results. 

The Conjugate ANOVA indicated that there was a 93% probabihty that the BNN 
approach had a smaller expected MSE than the Schnute approach. The Conjugate ANOVA 
estimated an expected MSE difference between the two techniques of 180 mm^. This trans- 
lates to approximately 1.7 m^ of timber per hectare. 

In the experiments conducted in this report, the factors that have been varied are 
the initial tree density and thinning strategy. Other relevant factors that have remained 
constant are the tree species, the forest plot location and sample ages that have been used 
for the training and test plots. These factors would need to be varied to ascertain which 
method was better in general. 

Another extension to this report, that might be of interest, is a comparison of Bayesian 
credibility intervals and Frequentist confidence intervals for NN modeling. It also might be 
of interest to compare how other 'nonparametric' methods, such as splines, would do at 
forest tree growth curve modeling. 

One of the broader issues addressed has been on the use of very complex techniques 
such as BNNs for relatively simple nonlinear regression problems. As has been demonstrated 
BNNs can compare favourably relative to more parametric, traditional approaches. The 
availability of free software removes, to some extent, implementation difficulties. However, 
it is my belief that for such complex techniques to gain more general acceptance, the choice 
of hyper-priors and evaluation of MCMC convergence needs to be made more automatic. 
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Appendix A 



Data 



The data used in this research report consisted of 12 plots of Pinus roxburghii Sargent, a pine 
native to the Himalayas. The plots were planted at an espacement of 1.8 by 1.8 m^ at the 
Weza forest plantation in what is now known as Kwa-Zulu Natal, South Africa. Each plot 



covered 0.08 ha and was surrounded by a 29 m wide buffer zone of trees ( Falkenhagen 1997 ). 
The plots were part of a series of correlated curve trend (CCT) experiments established by 



O'Connor in 1935 (Bredenkamp 1984). Different thinning regimes were used for each plot. 
The average diameter at breast height (DBH) measurements for each plot are shown 



in Table A.l. The number of trees within each plot is shown in Table A. 2 
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Ago 


Plot 


1 1 2 1 4 1 6 1 8 1 9 1 10 1 11 1 12 1 15 1 16 | 18 






5 


32.2 


39.3 


46.4 


40.2 


40.6 


27.8 


27.4 


33.1 


47.9 


43.9 


58.3 


49.0 


6 


52.1 


64.2 


74.1 


66.5 


68.1 


49.0 


47.8 


59.3 


75.7 


67.7 


87.8 


73.5 


8 


87.0 


109.0 


125.8 


119.3 


119.3 


85.6 


85.8 


104.1 


125.5 


107.7 


143.7 


116.9 


9 


99.3 


127.1 


152.5 


140.7 


143.7 


101.6 


109.8 


133.0 


145.3 


121.9 


163.0 


145.3 


10 


111.5 


142.7 


173.3 


166.6 


170.9 


116.5 


126.9 


154.4 


182.1 


135.2 


184.6 


166.2 


11 


120.3 


155.9 


190.1 


191.3 


195.1 


127.5 


140.3 


170.4 


202.4 


144.7 


201.9 


185.9 


12 


124.7 


161.7 


200.4 


200.8 


208.4 


130.3 


144.0 


177.6 


211.8 


147.5 


207.8 


192.9 


14 


135.1 


179.0 


226.4 


231.6 


273.9 


161.4 


179.2 


218.9 


258.3 


206.4 


244.1 


223.0 


16 


142.7 


194.1 


250.7 


262.1 


292.1 


176.7 


201.9 


248.9 


295.8 


236.2 


276.9 


250.3 


IS 


149.8 


205.0 


267.7 


286.5 


324.6 


206.5 


242.0 


281.9 


322.3 


257.7 


302.0 


265.7 


20 


157.2 


213.4 


284.5 


308.5 


354.1 


223.5 


262.7 


308.7 


347.6 


281.0 


320.6 


279.5 


22 


164.1 


221.4 


297.9 


327.2 


387.1 


254.5 


298.2 


335.8 


370.7 


299.7 


338.7 




24 


166.4 


227.6 


308.4 


341.7 


413.5 


272.0 


319.2 


357.1 


385.4 


314.9 


353.0 


305.6 


26 


172.0 


234.2 


317.9 


354.6 


437.1 


303.5 


338.3 


377.2 


400.9 


328.3 


366.0 


318.8 


28 


176.3 


240.4 


328.5 


367.5 


461.0 


325.1 


358.3 


397.5 


417.2 


342.0 


378.4 


333.4 


30 


180.4 


243.9 


334.8 


377.5 


477.3 


343.7 


373.0 


413.0 


425.4 


349.9 


385.2 


342.9 


32 


183.2 


248.6 


339.9 


384.9 


495.3 


358.4 


384.4 


428.4 


435.2 


357.7 


394.7 


353.2 


34 


189.2 


251.8 


348.6 


395.6 


513.4 


385.2 


385.4 


442.8 


441.7 


366.0 


405.6 


364.9 


37 


202.8 


256.2 


352.6 


401.7 


530.6 


403.2 


402.0 


459.2 


452.5 


374.7 


412.3 


377.3 



Table A.l: The average diameter at breast height (DBH) data. The age is measured in years 
and the DBH in mm. 



Ago 


Plot 


1 1 i 1 4 1 6 1 8 1 9 1 lo 1 11 1 li 1 IS 1 16 1 18 






5 


226 


120 


120 


120 


120 


233 


232 


160 


80 


214 


80 


80 


6 


226 


120 


120 


120 


120 


233 


232 


160 


80 


214 


80 


80 


8 


226 


120 


80 


80 


80 


233 


232 


160 


80 


214 


80 


80 


9 


226 


120 


60 


60 


60 


233 


162 


80 


40 


214 


80 


41 


10 


226 


120 


60 


40 


40 


233 


162 


80 


40 


213 


80 


41 


11 


226 


120 


60 


30 


30 


233 


162 


80 


40 


213 


80 


41 


12 


226 


120 


60 


30 


20 


233 


162 


80 


40 


213 


80 


41 


14 


226 


120 


60 


30 


20 


163 


80 


40 


20 


40 


39 


21 


16 


200 


120 


60 


30 


10 


162 


80 


40 


20 


40 


39 


21 


18 


200 


120 


60 


30 


10 


82 


39 


20 


20 


40 


39 


21 


20 


200 


120 


60 


30 


10 


80 


39 


20 


20 


40 


39 


21 


22 


200 


120 


60 


30 


10 


40 


20 


20 


20 


40 


39 


21 


24 


200 


120 


60 


30 


10 


40 


20 


20 


20 


40 


39 


21 


26 


200 


119 


60 


30 


10 


20 


20 


20 


20 


40 


39 


21 


28 


200 


119 


60 


30 


10 


20 


20 


20 


20 


40 


39 


21 


30 


200 


118 


59 


30 


10 


20 


20 


20 


20 


40 


39 


21 


32 


200 


117 


59 


30 


10 


20 


20 


20 


20 


40 


39 


21 


34 


193 


117 


59 


30 


8 


17 


19 


20 


19 


40 


39 


20 


37 


190 


117 


59 


30 


8 


17 


19 


20 


19 


40 


39 


20 



Table A. 2: The number of trees within a plot. The age is measured in years. 



70 



Bibliography 



J. Bergcr. Statistical Decision theory and Bayesian Analysis. Springer, 1985. 

J. M. Bernardo and A. F. M. Smith. Bayesian Theory. John Wiley & Sons, 1994. 

C. M. Bishop. Bayesian methods for neural networks. Technical Report NCR/95/009, De- 
partment of Computer Science and Applied Mathematics, Aston University, Birmingham, 
B4 7ET, U.K. http://www.ncrg.aston.ac.uk, 1995. 

G. E. Box and G. C. Tiao. Bayesian Inference in Statistical Analysis. John Wiley and Sons, 
Inc., 1992. 

B. V. Bredenkamp. The C.C.T. concept in spacing research. In Symposium on Site and 
Productivity of Fast Growing Plantations, pages 313-332, Pretoria and Pietermaritzburg, 
South Africa, april 1984. 

B. V. Bredenkamp and H. E. Burkhart. An examination of spacing indices for Eucalyptus 
grandis. J. For. Res, 20:1909-1916, 1990. 

B. V. Bredenkamp and T. G. Gregoire. A forestry application of Schnute's generalized 
growth function. Forest Science, 34(3):790-797, 1988. 

S. Brittain and L. M. Haines. Nonlinear models for neural networks. In Mathematics of 
Neural Networks: Models, Algorithms and Applications, pages 129-133. Kluwer, 1997. 

W. Buntine and A. Weigend. Bayesian back-propagation. Complex Systems, 5:603-643, 
1991. 

B. Cheng and D. M. Titterington. Neural networks: A review from a statistical perspective. 
Statistical Science, 9(l):2-54, 1994. 

M. K. Cowlcs and B. P. Carlin. Markov chain Monte Carlo convergence diagnostics: A com- 
parative review. Available from ftp://muskic.biostat.umn.edu/pub/1994/rr94-008.ps.Z, 
1995. 

G. Cybenko. Approximation by superpositions of a sigmoidal function. Mathematics of 
Control, Signals and Systems, 2:303-314, 1989. 

R. L. Czaplewski, R. M. Reich, and W. A. Bcchtold. Spatial autocorrelation in growth of 
undisturbed natural pine stands across Georgia. Forest Science, 40(2):314-328, May 1994. 

71 



G. A. Dover. rDNA world falling to pieces. Nature, 336:623-624, 1988. 

N. R. Draper and H. Smith. Applied regression analysis. John Wiley and Sons, New York, 
1981. 

S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth. Hybrid monte carlo. Physics 
Letters B, 195:216-222, 1987. 

M. Evans and T. Swartz. Methods for approximating integrals in statistics with special 
emphasis on Bayesian integration problems. Statistical Science, 10(3):254-272, 1995. 

E. R. Falkenhagen. Growth modeling of Pinus roxburghii in South Africa. Southern African 
Forestry Journal, (178):31-38, March 1997. 

A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin. Bayesian Data Analysis. Chapman 
and Hall, 1995. Lots of detailed examples, including MCMC and hierarchical modeling. 

S. Geman, E. Bienenstock, and R. Doursat. Neural networks and the bias/variance dilemma. 
Neural Computation, 4:1-58, 1992. 

W. R. Gilks, S. Richardson, and D. J. Spiegelhalter. Markov Chain Monte Carlo in Practice. 
Chapman & HaU, 1996. 

P. H. Goodman. NevPropS User Manual. University of Nevada, URL: 

ftp://unssun.scs.unr.edu/pub/goodman/nevpropdir/, 1996. 

C. Gordon. Predicting forest tree growth using a neural network. In C. J. Wright, editor. 
Proceedings of the Twentieth South African Symposium on Numerical Mathematics, pages 
28-45. SANUM and the Department of Computer Science, University of Natal, Durban, 
South Africa, Department of Computational and Applied Mathematics, University of the 
Witwatersrand, Johannesburg, South Africa, July 1994a. 

C. Gordon. The use of cross-validation in neural network extrapolation of forest tree growth. 
In Proceedings of the Pattern Recognition Association of South Africa, 1994b. 

E. J. Green, J. Francis A. Roesch, A. F. M. Smith, and W. E. Strawderman. Bayesian esti- 
mation for the three-parameter WeibuU distribution with tree diameter data. Biometrics, 
50:254-269, March 1994. 

E. J. Green and W. E. Strawderman. A comparison of hierarchical Bayes and empirical 
Bayes methods with a forestry application. Forest Science, 38(2):350-366, April 1992. 

E. J. Green and W. E. Strawderman. A Bayesian growth and yield model for slash pine 
plantations. Journal of Applied Statistics, 23(2):285-299, 1996. 

B. T. Guan and G. Z. Gertner. Machine learning and its possible roles in forest science. 
Artificial Intelligence Applications in Natural Resource Management Journal, 5(2):39-52, 
1991a. 



72 



B. T. Guan and G. Z. Gcrtncr. Modeling red pine tree survival with an artificial neural 
network. Forest Science, 37(5):1429-1440, 1991b. 

B. T. Guan and G. Z. Gertner. Using a parallel distributed processing system to model 
individual tree mortality. Forest Science, 37(3):871-885, 1991c. 

B. T. Guan and G. Z. Gertner. Modeling individual tree survival probability with a random 
optimization procedure: An artificial neural network approach. Artificial Intelligence 
Applications in Natural Resource Management Journal, 9(2):39-52, 1995. 

B. Hassibi and D. Stork. Second order derivatives for network pruning: Optimal brain 
surgeon. In C. Giles, S. Hanson, and J. Cowan, editors, Advances in Neural Information 
Processing Systems 5, pages 164-171, San Mateo, Galifornia, 1993. Morgan Kaufmann. 

S. Haykin. Neural Networks: A Comprehensive Foundation. Macmillan College Publishing 
Company, New York, 1994. 

K. Hornik, M. Stichcombe, and H. White. Universal approximation of an unknown mapping 
and its derivatives using multilayer feedforward networks. Neural Networks, 3:551-560, 
1990. 

E. T. Jaynes. Probability Theory: The Logic of Science. 1996. The latest version is available 
by ftp from bayes.wustl.edu. 

W. Jefferys and J. Berger. Ockham's razor and Bayesian analysis. American Scientist, 80: 
64-72, 1992. 

R. H. Jones and L. M. Ackerson. Serial correlation in unequally spaced longitudinal data. 
Biometnka, 77(44) :721-731, 1990. 

S. Kirpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. Science, 
220:671-680, 1983. 

R. L. Korol, K. S. Milner, and S. W. Running. Testing a mechanistic model for predicting 
stand and tree growth. Forest Science, 42(2):139-153, 1996. 

A. M. Kshirsagar. Multivariate Analysis. Marcel Dekker, inc. New York, 1976. 

S. Y. Kung. Digital Neural Networks. PTR Prentice HaU, 1993. 

P. Lambert. Modeling of nonlinear growth curve on series of correlated count data measured 
at unequally spaced times: A full likelihood based approach. Biometrics, 52:50-55, Mar. 
1996. 

Y. Le Cun, J. S. Denker, and S. A. SoUa. Optimal brain damage. In Advances in Neural 
Information Processing Systems, pages 598-605. Morgan Kauffmann, San Mateo, CA, 
1990. 

D. J. G. MacKay. Bayesian Methods for Adaptive Models. PhD thesis, Galifornia Insti- 
tute of Technology, Pasadena, Galifornia, 1992a. Can be obtained over the Internet at 
ftp://131.111.48.8/pub/mackay/. 

73 



D. J. C. MacKay. A practical Bayesian framework for backpropagation networks. Neural 
Computation, 4(3):448-472, 1992b. 

D. J. C. MacKay. Hyperparameters: Optimise, or integrate out? In G. Heidbreder, editor, 
Maximum Entropy and Bayesian Methods, Santa Barbara 1993, Dordrecht, 1994. Kluwer. 

M. E. McDill and R. L. Amateis. Fitting discrete-time dynamic models having any time 
interval. Forest Science, 39(3):499-519, 1991. 

W. T. Miller III, R. S. Sutton, and P. J. Werbos, editors. Neural Networks for Control. MIT 
Press, 1990. 

M. L. Minsky and S. A. Papart. Perceptrons. MIT Press, 1969. 

M. L. Minsky and S. A. Papart. Perceptrons. MIT Press, 2nd edition, 1990. 

R. M. Neal. Probabilistic inference using Markov chain Monte Carlo methods. Technical 
Report CRG-TR-93-1, Dept. of Computer Science, University of Toronto, 1993. 

R. M. Neal. Bayesian Learning for Neural Networks. Lecture Notes in Statistics No. 118. 
Springer- Verlag, New York, 1996. 

J. A. Nelder. The fitting of a generalization of the logistic curve. Biometrics, pages 89-110, 
Mar. 1961. 

J. A. Nelder. Note: An alternative form of a generalized logistic equation. Biometrics, pages 
614-616, Dec. 1962. 

A. J. O'Connor. Forest research with special reference to planting distances and thinning. 
British Empire Forestry Conference, South Africa, 1935. 

J. S. A. Oshu. Matrix model for tree population projection in a tropical rain forest of 
South- Western Nigeria. Ecological Modeling, 59:247-255, 1991. 

A. Papoulis. Probability & Statistics. Prentice Hall, 1990. 

V. V. Phansalkar and P. S. Sastry Analysis of the back-propogation algorithm with mo- 
mentum. IEEE Transactions on Neural Networks, 5:505-510, 1994. 

S. J. Press. Bayesian Statistics. Wiley, 1989. 

K. J. Puettmann. Evaluation of the size-density relationships for pure red Alder and 
Douglas-Fir stands. Forest Science, 39(l):7-27, Feb. 1993. 

C. E. Rasmussen. Evaluation of Caussian Processes and Other Methods for Non-linear 
Regression. PhD thesis. University of Toronto, 1996. 

D. A. Ratowsky. Nonlinear Regression Modeling. Marcel Dekker, 1983. 

M. R. Reynolds, T. E. Burk, and W.-C. Huang. Goodness-of-fit tests and model selection 
procedures for diameter distribution models. Forest Science, 34(2):373-399, 1988. 



74 



F. J. Richards. A flexible growth function for emperical use. Journal of Experimental 
Botany, 10:290-300, 1959. 

B. D. Ripley. Neural networks and related methods for classification. Journal of the Royal 
Statistical Society B, (3):409-456, 1994. 

D. E. Rumelhart, J. L. McClelland, and the PDP Research Group. Parallel Distributed 
Processing: Explorations in the Micro structure of Cognition. MIT Press, 1986. 

S. Saarinen, R. Braniley, and G. Cybenko. Ill-conditioning in neural network training prob- 
lems. SI AM Journal on Scientific Computing, 14(3):693-714, 1993. 

W. S. Sarlc. Stopped training and other remedies for overfitting. In Proceedings of the 27th 
Symposium on the Interface, 1995. 

J. Schnute. A versatile growth model with statistically stable parameters. Canadian Journal 
of Fish and Aquatic Science, 38:1128-1140, 1981. 

G. A. F. Seber and C. J. Wild. Nonlinear Regression. John Wiley & Sons, 1989. 

P. Soares, M. Tome, J. P. Skovsgaard, and J. K. Vanclay. Evaluating a growth model for for- 
est management using continuous forest inventory data. Forest Ecology and Management, 
71:251-265, 1995. 

D. Spiegelhalter, A. Thomas, N. Best, and W. Gilks. BUGS 0.5 Examples, Volume 1 (version 
\). MRC Biostatistics Unit, Institute of Public Health, Robinson Way, Cambridge CB2 
2SR, URL: http://www.mrc-bsu.cam.ac.uk/bugs/mainpage.html, 1996. 

H. H. Thodberg. Ace of Bayes: application of neural networks with pruning. Technical 
Report 1132 E, Danish meat research institute, 1993. 

J. K. Vanclay. Sustainable timber harvesting: Simulation studies in the tropical rainforests 
of North Queensland. Forest Ecology and Management, 69:299-320, 1994. 

J. K. Vanclay. Growth models for tropical forests: A synthesis of models and methods. 
Forest Science, 41(l):7-42, Feb. 1995. 

J. K. Vanclay, J. P. Skovsgaard, and C. P. Hansen. Assessing the quality of permanent 
sample plot databases for growth modeling in forest plantations. Forest Ecology and 
Management, 71:177-186, 1995. 

A. S. Weigend, D. E. Rumelhart, and B. A. Huberman. Generalization by weight-elimination 
with applications to forecasting. In R. P. L. et. al., editor. Advances in Neural Information 
Processing Systems 3, pages 875-882. Morgan Kaufmann, 1991. 

P. M. Williams. Bayesian regulerization and pruning using a Laplace prior. Neural Compu- 
tation, 7:117-143, 1995. 

R. L. Winkler. An Introduction to Bayesian Inference and Decision. Holt, Rinehart and 
Winston, inc., 1972. 



75 



B. Zcidc. Analysis of growth equations. Forest Science, 39(3):594-616, Aug. 1993. 

A. Zcllncr. The Bayesian method of moments (BMOM) theory and apphcations. Un- 
pubhshed report, HGB Alexander Research Foundation, Graduate School of Business, 
University of Chicago, 1996. 



76 



